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1  •  Research  Project  Overview 


This  report  describes  the  progress  and  results  of  the  University 
of  Southern  California  image  processing  research  study  for  the  period  of 
1  September  1973  to  28  February  1974.  The  image  processing  research 
study  has  b«en  subdivided  into  five  projects: 

Image  Coding  Projects 

Image  Restoration  and  Enhancement  Projects 
Image  Data  Extraction  Projects 
Image  Analysis  Projects 
Image  Processing  Support  Projects 

In  image  coding  the  orientation  of  the  research  is  toward  the  development 
of  digital  image  coding  systems  that  represent  monochrome  and  color  images 
with  a  nimmal  number  of  code  bits.  Image  restoration  is  the  task  of  im¬ 
proving  the  fidelity  of  an  image  in  the  sense  of  compensating  for  image  de  ¬ 
gradations.  In  image  enhancement,  picture  manipulation  processes  are 
performed  to  provide  a  more  subjectively  pleasing  image  or  to  convert  the 
image  to  a  iorm  more  amenable  to  human  or  machine  analysis.  The  objec¬ 
tives  of  the  image  data  extraction  projects  are  the  registration  of  images, 
detection  of  objects  within  pictures  and  measurements  of  image  features. 

Tl'«  image  analysis  projects  comprise  the  background  research  effort 
into  the  basic  structure  of  images  in  order  to  develop  meaningful  quantita¬ 
tive  characterizations  of  an  image.  Finally,  the  image  support  projects 
include  research  on  image  processing  computer  languages  and  the  develop¬ 
ment  of  experimental  equipment  for  the  sensing,  processing,  and  display 
of  images. 

The  next  section  of  this  report  summarizes  some  of  the  research 
project  activities  during  the  past  six  months.  Sections  3  to  7  describe  the 
research  effort  on  the  projects  listed  above  during  the  reporting  period. 
Section  8  is  a  list  of  publications  by  project  members. 


2.  Research  Project  Activities 


Significant  research  project  activities  of  the  past  six  months  are 
summarized  below: 

Alexander  A.  Sawchuk  has  been  appointed  editor  of  the  May- June 
1974  issue  of  Optical  Engineering,  the  journal  of  the  Society  of  Photo-Optical 
Instrumentation  Engineers.  This  is  a  special  issue  devoted  to  optical  and 
digital  image  and  information  processing,  and  will  contain  more  than  eleven 
papers  discussing  theory  and  applications  of  such  systems.  Topics  in  this 
issue  include:  hybrid  optical/digital  systems  using  real-time  input/output 
devices;  digital  holograms  for  optical  processing;  and  image  processing  of 
synthetic  aperture  data. 

Harry  C.  Andrews  has  been  appointed  guest  editor  of  the  May  1974 
issue,  of  Computer  .  the  journal  of  the  IEEE  Computer  Society.  The  issue, 
entitled  "Computer  Image  Processing,"  consists  of  seven  papers  spanning 
picture  coding,  image  restoration,  and  digital  image  processing  facilities 
at  Aerospace,  EG&G,  and  AEC  facilities.  The  issue  will  include  photo¬ 
graphic  examples  of  computer  processed  images  both  in  black  and  white  as 
well  as  in  color.  In  addition  to  normal  circulation,  the  issue  is  being  over¬ 
printed  for  additional  distribution  at  the  NCC  conference  in  Chicago. 
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3.  Image  Coding  Projects 

The  research  effort  in  image  coding  has  been  directed  toward  a  wide 
variety  of  applications.  Coding  systems  are  under  investigation  for:  monochrome 
and  color  imagery;  alow  scan  and  real  time  television;  and  information  pre¬ 
serving  and  controlled  fidelity  operation.  Results  of  this  research  study 
during  the  past  six  months  are  summarized  here  and  presented  in  detail 
in  subsequent  sections. 

In  the  first  report  an  interpolative  data  representation  is  utilized  to 
develop  three  image  coding  algorithms.  One  algorithm  is  based  upon  spatial 
domain  coding,  another  upon  transform  domain  coding,  and  the  third  is  a 
hybrid  coding  scheme.  Computational  requirements  for  the  algorithms  are 

specified  and  the  image  performance  is  evaluated  for  several  pictorial 
examples. 

The  next  report  describes  the  analysis  of  logarithmic  quantization 
scales  for  monochrome  image  quantization.  Quantization  errors  are 
evaluated  in  terms  of  a  moael  of  the  human  visual  process. 

Spline  functions,  which  are  a  special  class  ot  truncated  polynomials, 
are  known  to  be  quite  accurate  for  the  approximation  of  one  dimensional 
functions.  Their  use  in  image  approximation  for  purposes  of  bandwidth 
reduction  is  explored  in  the  next  report. 

In  the  following  two  reports,  the  concept  of  transform  domain  spectrum 
extrapolation  and  interpolation  for  image  coding  is  investigated.  With  these 
techniques,  transform  coefficient  quantization  error  can  be  reduced  by 

post-processing  at  the  coder  with  a  significant  impi  cement  in  image 
quality. 


The  last  report  considers  extensions  to  the  universal  coding  concept. 

In  particular,  a  rate  distortion  bound  is  established  for  coding  image  sources 
with  unknown  probabilities. 


3.  1  Image  Coding  via  Two  Dimensional  Interpolative  Representations 
Anil  K.  Jain 

For  finite  discrete  signals,  non-causal  "interpolative"  representations 
may  be  used  for  coding.  These  non-causal  representations  lead  to  three 
different  coding  algorithms  in  the  spatial,  hybrid  and  frequency  domains. 

Interpolative  Modeling  For  simplicity  in  presentation,  only  the  first 
order  stationary  Markov  signal  will  be  considered.  Let  {uj,  i  =  0,  1,  .  . . , 

N,  N+l  represent  such  a  signal  with  zero  mean  and  autocorrelation  given  by 

E[u.u.]  =  p^1"^  (1) 

Representation  of  a  sequence  {uj  denotes  a  relationship 

:£{u. }  =  v.  (2) 

i  i 

such  that  the  sequence  {u.}  can  be  reconstructed  from  v^.  For  example 

u..  -  pu.  ~  v.  (3) 

i  +  l  i  i 

is  the  Markov  representation  of  eq.  (1)  with 

P2  I  E[vZl  =  (1-p2)  (4) 

i  i  ' 

It  can  be  shown  that  the  linear  minimum  mean  square  representation  of 
eq.  (1)  is  given  by 


u.  -  — (u.  ,  ,  +u.  )  =  v. 

2  i+l  i-l  i 


i  =  1 ,  .  . . ,  N 


1  1+p' 


uo  ■  pui  =  vo 


(5a) 

(5b) 


UN+1  "  PUN  VN+1 


(5c) 


This  representation  is  such  that  E[v  ]  is  minimum  compared  to  all  other 
linear  representations  and  is  given  by 
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pi. 


E[v2]  £  f)2  =  1LEJ 

!■  «  ..  b, 

(1+P  ) 


(6) 


i'For  a  discussion  of  the  correlation  properties  of  and  generalization  of 

the  above  statement  to  non- stationary  case  see  [l].)  Also  observe  that 
2  2 

<  since  0  <p  <  1. 

For  a  two  dimensional  zero  mean  discrete  signal  u  with 

ij 


E[u..u.  ]  =  plnl+l 

LiJi+n,j+mJ 


m 


"horizontal"  and  "vertical"  representations  are  defined  as 

ev(l# j)  =  u..  -  a(u.  .  ,  +u.  .  ) 

h  U  M+1  i,j-l 


E(i*j)  =  u..  -  a(u  .+u.  .) 

v  iJ  i+l,J  i-l,  j 


where  a  =  p/(l+p  ).  Then  the  representation 


(7) 


(8a) 

(8b) 


ij  ij  2  ^Ui+i,  j  +  Ui-1,  j+  Ui,  j+1  +  Ui,  j-1* 

l  v 


is  such  that  E[e£+G2l  is  minimized,  and 


(9) 


This  is  in  contrast  with  the  Markov  representation  of  eq.  (7)  given  by 


(9a) 


u. .  =  pu. 
V) 


with 


••  i  •  +  Pu.  .  .  -  p“u.  .  .  ,  +  e.. 
l_  1  *  J  i-l, J-l  ij 


2.2 


E[e..]  =  (1-p^) 


(9b) 

(9c) 


Comparison  between  eqs.  (9a)  and  (9c)  shows  that  E[e  l<E[e  i  for 

ij  J  L  ij  J 

0  <  p  s  o.  786.  For  values  of  p  ~  1,  the  two  values  have  small  mean  square 
difference.  Also,  in  eq.  (9),  the  coefficient  'a1  is  relatively  insensitive  to 
changes  in  image  statistics  (p  parameter).  In  fact 
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(10) 


da  (1-P?‘)  /d^\ 

a  (Hp2)^P/ 

Therefore  in  the  vicinity  of  p  =  1,  small  changes  in  im^ge  statistics  will 
not  alter  the  performance  of  the  interpolative  representation. 

Coding  and  Reconstruction  Algorithms  It  should  be  recognized  that 
if  the  e  become  observables  in  eq.  (9),  then  this  equation  has  to  be  solved 
(for  reasons  of  stability)  as  a  boundary  value  problem  with  end  conditions 

u0  j>  ^+1  y  ui  o  and  ui  N+l  ^noWn*  ^or  simplicity,  it  will  be  assumed 
here  that  these  quantities  are  zero,i.e.,  the  picture  under  consideration  is 
imbedded  in  a  dark  background.  One  might  question  the  validity  of  this 
assumption  in  view  of  the  stationary  statistics  of  eq.  (7);  however,  this 
assumption  is  actually  non-essential  for  the  coding  algorithms  below. 
Equation  (9)  can  then  be  rewritten  in  vector  form  as 

2  Uj+1  ‘  Quj  +  2  Uj-1  =  'Gj  ’  (H) 

where  u.  and  e.  now  represent  N  x  1  column  vectors  of  elements  (u  ,  .  . . ,  u 

J  J  lj  Nj 

and  (e ij,  . . . ,  respectively.  The  matrix  O  is  a  symmetric  tridiagonal 
matrix  given  by  the  elements 


[  1 

i  =  j 

Q..  =  J  -a/2 

ij  / 

lHl-i 

(12) 

\  0 

otherwise 

Algorithm  Al:  (Spatial  Domain  Coding,  Figure  la) 

1.  Quantize  and  code  e^,  after  obtaining  it  through  eq.  (9).  Let 

e'1  .  denote  the  received  signal. 

LJ 

2.  The  reconstructed  sequence  u’1'  is  obtained  by  solving 


V 

=  R.ur 
j  j 

+  s. 

J 

uo  *  0 

(13a) 

v. 

rn  =  0 

(13b) 

V 

-  Vi 

ft  i +  sj) 

CO 

2 

ii 

o 

(13c) 

-6- 


ORIGINAL 
IMAGE 


INTERPOLATE 

LI- 

,Jr/ 

MODEL 

_V 

QUANTIZER 


€M =  uii~g.i  =  uii-T  ( WuMj+,WuiJ  J 


* 

€ij 


RCV 


DECODER 


* 

U.. 


Jmodel 

a)  Spatial  domain  coding  algorithm  A1 


u 


ij 


i_ 

INTERPOLATIVE 

Uii.( 

MODEL 

€ . . 


A 

€  -  =T€. 
J  J 


QUANTIZER 


A.  * 

€«j 

•'N* 

UM 

u. 

'J 

DECODER 

■a 

UpTUj 

RCV 

JMODEL 

b)  Hybrid  domain  coding  algorithm  A2 


Figure  3. 1,1  Image  coding  algorithms. 


Eq.  (13)  can  be  solved  in  about  (N  log£  N)  computations, 
employing  the  structure  of  matrix  Q.  [l,  2,  3] 


Algorithm  A2;  (Hybrid  Domain  Coding,  Figure  lb) 

1.  Find  e.  =Te.,  where  T  is  an  N  X  N  matrix  containing  /~~~~ '  sin ) 
J  J  6  V  N+l  \N+1/ 

terms.  It  can  be  shown  [3"]  that  T  is  idempotent  and  diagonalizes 


Q  i.  e. 


TQT  =  A  = 


2.  Quantize  e  with  the  number  of  bits  allocated  for  the  i  row 

being  proportional  to  logfl  -  cos  |  \  to  obtain  e*  . 

\  N+l/  ij 

3.  The  reconstructed  image  u  is  obtained  bv 

ij 


i'-  ^  A  j’s 

u:  =  Tu. 


uij+l  =  rijuij +  8ij 
r..  1  =  7  (x-  -Tr  ) 

ij-1  2  \  i  2  ij/ 

j  =  r.  .  e*  +  s  \ 
V)-1  i » J  -  !\a  ij 


j  ,  A  .  ITT 

ana  =  1  -  a  cos  —  _ 

l  N+l 

2 

Again  eq.  (15)  requires  about  (N  log£  N)  computations. 

Algorithm  A3:  (Transform  Domain  Coding,  Figure  lc) 

1.  If  e  is  N  XN  matrix  of  element  e..,  then  first  find 

ij) 

A 

e  =  TeT 

A 

2.  Quantize  e..  by  allocai'ng  n  .  bits  to  it  such  that 

lJ  VI 


n..a  log(|a.  +  |i  ) 


where  U.  =  1  -  2a  c 

i 


°ste) 


i  =  1, 


(17) 

(17a) 
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3.  The  reconstructed  image  is  given  by 


lur.j  -  U'  =  TU  T 
lJ 

a  2e'\ 

and  u  =  - ii- 

ij  M.+M. 

i  J 

References  [l,  2"|  provide  details  and  generalizations. 


(18) 

(18a) 


Examples,  Implementation  and  Computational  Considerations  Algorithms 

A2  and  A3  have  been  simulated  for  the  256  x  256  pixel  girl  image  of  fig.  2a. 

The  average  value  of  p  for  this  picture  is  0.  96.  In  all  the  simulations  the 

actual  value  used  was  p  =  1.  The  difference  in  e  .  for  these  two  values 

2  lJ 

visually  and  quantitatively  (in  terms  of  £e..)  both  was  insignificant. 

Figure  2b  shows  the  display  cf  |e  |.  Figure  2c  contains  the  encoded  image 
according  to  algorithm  A2  with  3  bits /pixel  on  the  average  using  a  uniform 
quantizer.  The  entropy  of  the  quantized  signal  §*  was  actually  2.  35  bits,  so 
that  a  variable  length  Huffman  code  could  be  employed  to  obtain  the  same 
image  with  2.35  bits/pixel.  Figure  2d  shows  the  quantized  |e*  |  signal 
according  to  algorithm  A3  and  bit  rate  of  1  bit/pixel.  Figures  2e  and  2f 
show  the  encoded  images  for  1  and  1.47  bits  /pixel. 


It  can  be  shown  that  the  total  computational  load  in  each  algorithm  is 
of  the  same  order.  In  algorithm  Al,  the  memory  and  computational  require¬ 
ments  on  the  transmitter  are  minimal  and  the  major  computational  burden  is 
at  the  receiver.  In  algorithm  A2,  the  total  processing  burden  is  roughly 
divided  in  a  1:2  ratio  (the  transmitter  needs  to  take  a  one  sided  transform 
and  the  receiver  solves  scalar  interpolative  equations  and  takes  an  inverse 
transform).  In  algorithm  A3,  the  processing  load  is  roughly  equally  divided 
between  the  transmitter  and  the  receiver.  Thus  the  three  algorithms  obtained 
via  a  single  representation  spell  out  three  different  communication  system 
architectures,  and  their  relative  use  therefore  depends  on  the  particular 
application.  In  this  sense  the  interpolative  representation  leads  to  a 
unification  of  some  of  the  different  methods  of  image  coding  employed 
currently. 
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(c)  hybrid  domain  encoded  image 


A 

(d)  quantized  signal  |  G  *  | 


(f)  encoded  image  1.47  bits /pixel 


(e)  encoded  image  1  bit/pixel 


Figure  3.  1-2.  Image  coding  results  via  interpolative  model. 


Finally,  it  should  be  remarked  that  the  representation  used  here  is 
but  one  member  of  a  class  of  similar  non-causal  and  semicausal  repre¬ 
sentations.  The  representation  reported  here  corresponds  to  a  discrete 
version  of  the  Poisson  equation  (\?  u=e).  The  relative  merits  of  other 
representations  is  currently  under  study  and  will  be  reported  in  the  future. 
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3.2  Optimal  Logarithmic  Quantization  for  Picture  Processing 
Francis  Kretz  and  Werner  Frei 

An  optimal  quantization  law  for  image  intensities  for  television  monitor 
display  of  digitized  and  prc^ssed  images  has  been  considered.  AIbo,  the 
effects  of  "brightness"  adjustments  for  televisiondiBplays  has  been  analyzed. 

Subjective  Criterion  A  possible  subjective  criterion  for  quantization 
distortion  is  to  postulate  that  the  decision  and  reconstruction  levels  should 
be  perceptually  equi-distant.  Since  it  iB  well  known  that  the  perception  of 
intensity  is  a  concave, monotomically  increasing  function  of  light  intensity 
A(I),  an  optimal  quantizer  in  the  above  sense  can  be  derived  from  A(I)  as 

Q(I)  =  int[(2N- 1  )A  (I/I  )]  (1) 

where  int[*  ]  denotes  the  nearest  integer  (OsQs  2N-1)  of  the  argument, 

N  is  the  number  of  bits  of  the  quantizer,  and  A(I)  is  normalized  so  that 

A(0)  =  0,  A(l)  =  1 
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Various  functions  for  A (I)  such  as  the  square  and  cube  root,  polynomials 
and  logarithms  have  been  proposed  to  fit  experimental  perceptual  data  [l]. 
The  parameters  and  range  of  validity  of  these  functions  depend  very  much 
upon  the  experimental  conditions  under  which  the  data  was  obtained.  In 
particular,  background  illumination  has  a  strong  influence  on  A(I). 

In  the  case  of  t  Revision  displays,  the  range  of  intensities  is  fairly 
well  defined  (about  two  orders  of  magnitude),  but  the  background  illumination 
for  each  pixel  is  a  complex,  more  or  less  random  field  (the  image  itself). 
Since  one  desires  to  design  one  quantizer  for  all  pixels  of  an  image,  the 
parameters  of  A(I)  will  be  the  result  of  some  compromise. 

Experiments  have  been  carried  out  to  determine  the  optimal  slope 
at  the  origin  of  the  A-function 


A(I)  =  b  loglQ(l  +  I/a)  (2) 

where 

b  =  (logio(1  +  a*)_1 

with  A(0)  =  0  and  A ( 1 )  =  1  .  For  this  A-function  the  slope  at  the  origin  is 


A'(0)  = 


a  loge(lO) 


In  the  first  set  of  experiment,  a  ramp  of  intensities  was  quantized 
according  to  eqs.  (1)  and  (2),  generating  a  set  of  grey  scales  with  A 1  (0 )  = 

4.4,  6.3,  8.0,  11.3,  16.1,  19.7,  22.9,  29.  9  (from  top  to  bottom,  respectively 
in  figure  la  (N  =  4)  and  figure  lb  (N  =  3  and  5  bits).  Then  the  parameter 
A'(0)  was  chosen  corresponding  to  the  scale  with  the  most  uniform  subjective 
spacing  of  intensities  over  the  entire  range.  It  should  be  pointed  out  that  the 
figures  reproduced  here  have  been  subjected  tc  a  series  of  distortions 
inherent  to  the  lithographic  process.  Several  observers  viewing  the  TV 
monitor  preferred  A'(0)  =  16  (slightly  steeper  than  the  Richter  scale,  see  [ll). 
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(a)  16  -  level  logarithmic  grey  scale 


In  the  second  experiment,  two  images  (the  SMPTE  "couple"  and  "girl") 
were  coarsely  quantized  (N  =  4  bits  per  pixel)  in  order  to  verify  both  the 
postulate  regarding  the  quality  criterion  and  the  optimal  slope  determined 
in  the  previous  experiment.  Figures  2  and  3  show  the  respective  original 
images  and  a  set  of  logarithmically  quantized  versions  with  different  para¬ 
meters  A'(0).  For  comparison,  linearly  quantized  versions  with  four  and 
five  bits  are  included. 

The  second  experiment  shows  that  the  optimal  parameter  A'(0)  depends 
on  the  picture  content.  Note  that  a  small  value  of  A'(0)  tends  to  create  large 
subjective  increments  in  dark  areas  and  vice-versa.  It  is  observed  that  the 
very  small  optimal  value  of  A 1  (0 )  for  the  "girl"  picture  is  a  consequence  of 
the  unusually  rare  occurrence  of  low  intensities.  In  fact,  the  histogram  of 
the  "girl"  picture  has  a  maximum  at  about  20%  intensity.  The  "couple" 
picture  has  a  more  typical  negative  exponential-like  histogram  and  the 
optimal  slope  A'(0)  =  16  is  the  same  as  determined  in  the  first  experiment. 

Minimum  square  error  quantization  laws  for  sources  with  given  ampli¬ 
tude  probabilities  have  been  studied  [2,  3*1.  Assuming  a  probability  density 
function  of  the  form  p(I)  «  k  exp(-al),  the  MSE  quantization  law  follows  a 
concave  monotomically  increasing  function  of  I  quite  similar  to  A (I).  The 
agreement  between  both  criteria  is  presently  being  studied  in  more  detail. 

Practical  Considerations  In  typical  computer  imag.i  processing 
environments,  images  are  usually  scanned  and  quantized  linearly.  It 
obviously  makes  little  sense  to  coarsely  requantize.  The  results  of  this 
study  are  therefore  primarily  relevant  to  scanner  or  coding  hardware. 

Comparing  the  quality  of  linearly  and  logarithmically  quantized  images, 
one  sees  that  one  bit  at  least  out  of  five,  possibly  two  out  of  eight  can  be  saved 
with  appropriate  quantization,  which  represents  a  20-25%  bandwidth  reduction, 
or  storage  saving,  whichever  is  relevant. 
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(b)  PCM  5  bits  linear 


(c)  PCM  4  bits  linear 


(d)  PCM  log.  4  bits  A'(0)=6.  3 


(e)  PCM  log.  4  bits  A'(0)=16.  1  (f)  PCM  log.  4  bits  A'(0)=22.  9 


Figure  3.2-2.  Examples  of  grey  scale  quantization  of  "couple"  picture 


L  a  fi 

\  v  I  m 

k  1 

In  the  context  of  hardware  considerations,  the  influence  of  the  "bright¬ 
ness"  control  of  a  television  monitor  has  been  analyzed.  Let  u  be  the  signal 
applied  to  the  grid  of  the  cathode  ray  tube  (CRT),  UQ  the  CRT  cut-off  voltage, 
Uj  tbe  bias  voltage  (controlled  by  the  "brightness"  setting)  and  v  the  intensity 
to  be  displayed.  The  light  intensity  of  the  CRT  is  proportional  to 

I  =  (u  +  UQ  +  u+U^+Uj  >0 

I  =  0  u+Uq+Uj  so 

The  Y  correction  is  carried  out  by  letting 

« •  v1/Y 

where  is  a  constant.  An  incorrect  setting  of  the  brightness  control 
VU1  =  ^  ^  0  does  not  simply  reduce  the  useful  dynamic  range  of  the 
display.  It  also  upsets  the  linearity  of  the  gamma-corrected  transfer 
function 

I  =  (kjV1/Y  +  AU)Y 

oince  the  eye  is  very  sensitive  to  errors  at  low  intensities,  the  effect  of 
A  J  is  quite  severe  for  the  lower  levels.  Figure  4a  shows  the  overall  transfer 
ft  .icvion  A(w)  where  w  represents  the  physiologically  companded  intensities 

v  =  a  [exp(^-)  -  1] 

and 

A  =  b  log(l  +  I/a) 

as  shown  in  figure  4b.  These  results  indicate  the  importance  of  brightness 
adjustments.  A  well  designed  hardware  system  should  have  an  automatic 
video  clamping  circuit  controlled  by  the  D/A  converter,  such  that  AU  is  set 
to  zero  when  a  digital  zero  is  read. 
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(a)  influence  of  brightness  control  on  overall  perceptual  system. 


(b)  perceptual  system  block  diagram. 


Figure  3.2-4.  Subjective  effect  of  brightness  control. 
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3.3  Image  Data  Compression  Using  Spline  Function! 
Faramaz  Davarian 


Spi...e  functions  can  be  used  to  simultaneously  compress  and  interpolate 
a  given  set  of  data.  Among  different  sets  of  bases  spanning  the  space  of 
splines,  the  B-spli„es  are  most  suitable,  since  they  possess  a  local  basis 
property  and  result  in  matrices  which  are  easily  invertable. 

A  dimensionality  and  subsequent  bandwidth  reduction  can  be  achieved  by 
a  least  squares  fit  of  n  points  using  m  (m  =-n)  basis  functions.  In  essence, 
the  data  compression  method  is  simply  a  transformation  of  the  n-dimensiolal 
space  of  data  points  into  a  smaller  m-dimensional  spline  space.  It  is  note¬ 
worthy  that  elements  of  the  m-dimensional  spline  space  will  directly  generate 
the  continuous  estimate  of  the  original  signal  rather  than  the  sampled  estimate. 

The  method  is  described  below  along  with  a  study  of  the  statistical  properties 
of  the  transform  domain. 


Definition 
of  real  numbers 
knots  tj,t  ,  .  .  .  , 
two  pioperties: 


_of_ Spline  Functions  Given  a  strictly  increasing  sequence 
,  tj ,  t2, . . . ,  tn,  a  spline  function  S(x)  of  degree  m  with 
tn  is  a  function  defined  on  the  real  line  having  the  following 


i) 

ii) 


In  each  interval  (t.f  t  )  for  i  =  0,1,... 
polynomial  of  degree  m  or  less. 

S(x)  an.!  its  derivatives  of  order  1,2,.. 


,  n-1,  S(x)  is  given  by 


.  ,  m- 1  are  continuous 


some 


on 


Thus,  a  spline  function  is  a  piecewise  polynomial  function  satisfying  certain 
conditions  regarding  continuity  of  the  function  and  its  derivatives. 
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It  is  generally  believed  that  in  many  circumstances  a  spline  function 
is  a  more  adaptable  approximating  function  than  a  polynomial.  This  is  based 
in  part  on  actual  numerical  experience,  and  in  part  on  mathematical  demon¬ 
strations  that  solutions  of  a  variety  of  problems  of  best  approximation  turn 
out  to  be  spline  functions.  A  spline  function  may  be  defined  in  terms  of  a 
truncated  power  function 


:-*V  -  { 


(x-t) 


(X  >t) 


(X  St) 


It  is  easily  seen  that  [l  ]  any  spline  of  degree  1  with  knots  t.,t  ,  .  »t 

1  2  n 

has  a  unique  representation  of  the  form 


S(x)  =  P(x) 


+  §  ci<x  -  V 


where  P(x)  is  a  polynomial  of  degree  l  or  less.  The  above  representation 
of  a  spline  function  normally  results  in  an  ill  conditioned  set  of  line.'  r 
equations.  To  overcome  this  deficiency  it  is  possible  to  introduce  a  new 
set  of  local  basis  functions  for  the  space  of  splines  [2"]. 

The  preceeding  considerations  lead  to  the  B- splines  studied  by 
Schoenberg  [3*1,  which  are  in  a  sense,  the  splines  of  minimal  support  for 
a  given  degree  (consisting  of  the  smallest  possible  number  of  intervals 
between  knots).  Figure  1  illustrates  the  typical  shape  of  a  B- spline  function. 

Data  Compression  by  Least  Squares  Method  Given  a  set  of  data  pairs 

(t.,y.)  for  i  =  1, 2,  .  .  .  ,  n,  which  can  be  interpolated  as  digitized  values  of  the 

points  of  the  graph  y  =  f(t),  let  the  unknown  function  f(x)  be  approximated  by 

a  linear  combination  of  suitably  chosen  functions  M.  (t),  M_(t), . . ,  ,  M  (t) 

12  m 

which  are  the  basis  splines. Then 


f(t)  =  c.  M  (t)  +  c  M  (t)  +  . . .  +c  M  (t) 
li  c  c.  mm 


M2, ...» Mm  form  a  complete  set  of  bases  for  the  space  of  m  data  points 
(m  dimensional).  This  basis  can  interpolate  m  data  elements  exactly. 


-20- 


Figure  3.3-1.  Typical  shape  of  a  B-spline. 


where  the  unknown  coefficients  Cj,  c^,  . . . ,  c^  are  independent  parameters 
to  be  determined,  and  m  ^n.  To  minimize  the  mean  square  approximation 
error 


Q  =  £(F(t.)-y.l2  =  [X£  c.M.(t.)  -  y.'l2 
1  1  1  L  i  j  J  J  t  ij 

differentiate  Q  with  respect  to  a^  and  set  the  result  to  zero  to  obtain 


3Q 
9  c. 


or 


=  -  y.}Mk(t.)  =  0 


£c  £  M .(t.)M,  (t.)  -  £y.M,  (t.)  =  0 
J  J  i  j  i  k  i  i  7i  k  i 


k  =  1,  2,  . . . ,  m 


In  matrix  form 


or 


where 


f?  MiMJ  K1  ■  [?  ^„'v] 


T  T 

B  BC  =  B  x 


Z  =  [y,,y2.y,....ynl 


£  *  fc,. 


c  1 
m  ’ 


B 


Mj(tj) 

M2(ti) 


M  (t.) 
m  1 


. M,(t_) 


1 


n 


Mz(V . 


M  (t  ) 
m  n 


with  M.(t  )  =  0,  for  |j  |  ^p.  Thus,  13  has  many  zero  value  off  diagonal 
entries.  The  vector  of  weighting  coefficients  is  then 
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rtfla 


iW 


where 


is  the  pseudoinverse  matrix  of  B.  Note  that  n  elements  of  data  vector  z 

are  mapped  into  m  elements  of  the  C  vector,  which,  represents  the  coefficients 
of  m  spline  basis  functions. 

The  estimate  of  f(t)  is  then 

A 

f(t)  =  c  M  (t)  +  c  M  (t)  +  ...+c  M  (t) 

c  l  mm 

A 

f(t.)  =  £  c.M.(t.)  i  =  1  n 

i  y  J  j'  1,  .  . .  ,n 

Let 

£  =  [F(tJ...F(t  )1T 

i  n 

£  =  BC  =  B^BfVy 
The  error  vector  e  can  then  be  expressed  as 

1  =  ft-  [±nX„-^BTB)-Yr]I 

Statistical  Properties  of  B-Spline  Coefficients  If  the  data  vector 
is  modeled  as  a  sample  of  a  vector  random  process  with  known  mean,  E{y), 
and  known  covariance,  K^,  the  B-spline  coefficients  given  by 

C  = 

are  also  random.  Their  mean  and  covariance  are 


E[C}  =  E  { =  B+E{y_) 

and 

E{CCT}  =  B+K  (B+)T 

—  - y  — 

If  the  data  vector  £  is  considered  a  sample  of  Markov  process  with  a 
correlation  coefficient  of  0  (0  <p  <  1)  between  each  adjacent  pixels  and 
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5 


10  15  16 

ELEMENT  NUMBER 

b)  N  =  16,  ELEMENT  CORRELATION  =.96 


a)  N  =  16,  ELEMENT  CORRELATION  =.95 


Figure  3.3-2.  Variance  of  B-spline  coefficients. 


self  correlation  coefficient  of  unity,  then 


K 


— y 


p  p 
i  p 


N- 1 
P 

N-2 

P 


1 


Figures  2a  and  2b  contain  two  plots  of  the  variance  function  of  c^  as  a  function 

of  i,  where  c.  is  the  i*^1  element  of  C.  The  plots  are  obtained  with  N  =  M  =  16. 
i  — 
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3.4  Positive  Extrapolation  of  Signals  and  Images 
Ali  Habibi,  Firouz  Naderi 

In  a  transform  coding  system  a  bandwidth  reduction  is  achieved  by 
discarding  a  number  of  transform  coefficients  of  a  natural  image.  Those 
coefficients  possessing  a  small  variance  are  of  low  information  content,  and 
replacing  them  by  zeros  at  the  receiver  results  in  a  rather  small  degradation 
in  the  quality  of  the  encoded  signal.  The  customary  approach  in  designing 
transform  coding  systems  has  been  to  substitute  the  missing  coefficients  at 
the  receiver  with  zeros.  However,  the  quality  of  the  coded  signal  improves 
by  extrapolation  of  the  missing  coefficients  from  those  which  have  been  trans¬ 
mitted.  The  problem  is  analogous  to  one  encountered  in  spectral  estimation 


of  data  where  the  covariance  function  is  first  estimated  for  a  number  of  lag 
values,  then  the  covariance  function  at  the  available  lag  values  is  extrapolated 
for  additional  lag  values  prior  to  taking  its  Fourier  transform.  The 
extrapolation  problem  as  applied  to  transform  coding  is  more  complicated. 

In  transform  coding  only  the  quantized  values  of  the  coefficients  are  available 
at  the  receiver;  thus,  one  is  forced  to  estimate  the  missing  components  of 
the  transformed  data  from  the  available  possibly  noisy  components. 

This  problem  has  been  analyzed  using  two  different  approaches.  The 
first  approach  is  statistical  and  is  based  upon  the  correlation  among  the 
transformed  components  using  suboptimal  transforms  such  as  the  Fourier, 
Hadamard,  and  Slant  transforms  .  The  second  approach  is  called 
positive  extrapolation  since  it  is  based  upon  the  positiveness  of  the  video 
data  and  the  fact  that  the  Toeplitz  matrix  constructed  from  the  Fourier 
coefficients  of  a  positive,  real  signal  is  always  positive  definite. 

Extrapolation  of  Signals  Let  F(0),  F(l),  . .  .  ,  F(N)  refer  to  the  first  N+l 
components  of  vector  F\  the  Fourier  transform  of  one  line  of  a  video  data, 
which  is  composed  of  M  points.  To  make  F^  real  one  must  generate  an  even 
function  by  first  reflecting  the  video  data  about  the  t  =  0  axis  and  then  taking 
a  Fourier  transform  of  the  even  signal. 

Now  consider  the  Toeplitz  matrix  T(N+1)  defined  as 

F(1 )  ...  F(N)  F(N+1 ) 

F(1 )  ...  F(N-l)  F(N) 

•  •  • 

•  •  • 

•  •  • 

F(N)  ...  F(1 )  F(0) 

Since  the  modified  video  data  is  positive  and  even,  and  the  real  matrix 
T(N+1)  is  semipositive  definite,  it  follows  that  the  determinant  of  T(N+1) 
as  a  function  of  F(N+1)  has  a  single  maximum.  Hence  the  allowable  values 
of  F(N+1)  are  those  that  make  the  determinant  of  T(N+1)  equal  to  zero  and 


\ 

/ 


T(N+1 )  = 


F(N+1 ) 
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all  values  in  between.  Expanding  the  determinant  of  T(n+1)  in  terms  of 
the  last  row  and  the  last  column  gives  an  expression  for  the  determinant  of 
T(N+1)  in  terms  of  F(0),  F(l),  . . . ,  F(N+1).  Since  F(0),  F(l),  . . . ,  F(N)  are 
known,  this  is  an  expression  for  the  determinant  of  ^(N+l)  in  terms  of 
F(N+1).  Choosing  F(N+1)  to  maximize  the  determinant  of  T(N+1)  gives  a 
recursive  algorithm  to  estimate  F(N+1)  from  F(0),  F(l),  . .  . ,  F(N).  The 
recursive  algorithm  can  be  used  further  to  estimate  F(N+2)  from  F(l),  .  .  . , 
F(N)  and  the  estimated  value  of  F(N+1)  i.e. 

N 

F(j)  =  A(k)F(  1  j-k  | )  for  j  =  N+l,...M-l  (2) 

k=l 

where  A(k),  k  =  1,  . .  . ,  N,  are  a  set  of  fixed  constants  specified  by  matrix 

T(N+1). 


Extrapolation  of  Images  The  positive  extrapolation  technique  discussed 
for  one-dimensional  signals  in  the  previous  section  can  be  generalized  to 
extrapolate  two-dimensional  spectral  density  functions  as  well  as  two- 
dimensional  Fourier  transform  of  images.  This  is  achieved  by  extending 
eq.  (2)  to  functions  of  two  variables  by  letting 


N  N 

F(i.j)  =  A(M)F(M|,  |j-k|) 


for  i,  j  =  0,  1,  . . . ,  M-l  (3) 


where  F(i,  j)  is  the  two-dimensional  discrete  Fourier  transform  of  the  image 

2  2 
and  consists  of  M  elements.  At  the  receiver  site  (N+l)  elements  are 

2 

available  and  these  (N+l)  elements  are  used  to  extrapolate  the  missing 
elements  prior  to  taking  the  inverse  Fourier  transform  to  obtain  a  re¬ 
construction  of  the  original  image.  Analogous  to  the  one-dimensional 
system,  the  original  picture  is  first  folded  along  the  x  =  0  and  y  =  0  axes 
to  generate  an  even  two-dimensional  array.  This  is  required  to  make  F(i,  j) 

an  array  of  real  elements.  Solving  eq.  (3)  for  A(k,<0)  is  straightforward 
2 

since  (N+l)  values  of  F(i,j)  are  known. 
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Experimental  Results  The  performance  of  the  positive  extrapolation 
methods  described  above  has  been  considered  for  a  number  of  examples. 

The  one-dimensional  example  is  a  discrete  signal  of  32  samples  that  consists 
of  a  pulse  superimposed  over  a  slowly  varying  background.  The  discrete 
Fourier  transform  of  this  signal  is  calculated  and  all  but  the  first  eight 
samples  substituted  by  zeros.  The  inverse  Fourier  transform  of  the  truncated 
signal  and  the  original  signal  are  shown  on  Figure  1.  The  truncated  signal 
in  the  transform  domain  is  then  extrapolated  to  recover  all  ?2  components. 

The  inverse  transform  of  the  extrapolated  signal  is  also  shown  on  Figure  1. 
The  extrapolated  signal  remains  positive  as  expected  and  tracks  the  original 
signal  significantly  closer  than  the  unextrapolated  signal. 

Figures  2  and  3  contain  examples  of  posit've  extrapolation  for  a  discrete 
(sin  x/x)(sin  y/y)  signal  and  a  block  pulse  signa. ,  respectively.  In  each  case 
the  Fourier  transform  has  been  truncated  preserving  the  transform  samples 
in  a  3  x  3  low  frequency  block  out  of  a  total  of  15  X  15  coetfi^ients.  The 
positive  extrapolation  process  is  seen  to  provide  a  significant  improvement 
over  the  reconstruction  without  extrapolation. 

3.  5  Transform  Domain  Spectrum  Interpolation 

Michael  N.  Huhns 

Quantization  occurs  whenever  continuous  physical  properties  are 
represented  numerically.  A  quantizer  is  a  zero-memory  nonlinear  device 
which  restricts  an  input  variable  to  a  fiiite  number  of  possible  output 
regions.  This  process  is  irreversible  and  information  is  invariably 
destroyed  since  only  the  region  containing  the  input  is  known  at  the  output. 
However  this  output  data  can  be  combined  with  a  priori  knowledge,  about 
the  input  to  reduce  the  amount  of  information  lost  by  interpolating  between 
the  discrete  outputs. 

In  transform  image  coding  a  block  of  image  pixels  undergoes  a  two 
dimensional  transformation  using  a  unitary  transform  such  as  the  Fourier, 


I 
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C  C  •  C  h T 


Positive  extrapolation  e: 


(a)  Original 


(c)  Extrapolated 

-3  .  The  original,  un- extrapolated  and  the  extrapolated 
two-dimensional  Signals.  Reduction  ratio  is-^  . 


Hadamard,  or  Slant  transform.  Next,  the  transform  coefficients  are 
quantized  and  coded  for  transmission.  Figure  1  illustrates  a  typical  bit 
assignment  for  a  zonal  quantization  and  coding  algorithm.  The  number  of 
quantization  levels  assigned  to  the  coefficient  at  coordinate  (u,  v)  is 


M(u,  v) 


2^(u»  v) 


(1) 


where  b(u,  v)  denotes  the  bit  assignment.  At  the  receiver,  the  quantized 
coefficients  are  reconstructed  and  an  inverse  transformation  is  performed 
to  obtain  an  image  estimate. 

If  a  transform  coefficient  is  quantized  to  zero  bits,  then  its  restoration 
is  equivalent  to  a  spectrum  extrapolation  as  outlined  by  Pratt  [l]  .  Those 
coefficients  that  are  quantized  to  two  or  more  levels  can  also  be  restored  by 
a  technique  called  spectrum  interpolation. 


Analysis  Le:  the  N  element  column  vector  x  with  probability  density 
p^(x)  denote  a  vector  of  input  data  samples.  For  two-dimensional  data  arrays, 
x  is  formed  by  column  scanning  the  data  array.  Each  data  sample  is  quantized 
into  one  of  M  output  regions,  denoted  by  D.,  i  =  1,2,..,,  M.  The  estimated 
value  of  xbased  upon  the  observed  D.  regions  is  the  quantizer  output  vector 
y..  The  average  error  in  this  estimate  is  then  defined  as 


S  = 


(2) 


where  e(*  )  is  an  arbitrary  error  weighting  criterion.  The  vector  of  estimates 
should  be  chosen  to  minimize  the  average  error.  This  choice  can  be 
determined  by  utilizing  the  principles  of  calculus  to  find  the  stationary 
points  of  the  error  surface  S  with  respect  to  each^.  Herce 


-J  [e(x-^)]p  (x)dx 

D.  -H 

i 


i=l  2, . .  . ,  M  (•)) 
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Figure  3.  5 
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1.  Typical  transform  domain  quantizing  bit  assignment. 
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with  the  assumption  that  the  error  function  e(*  )  is  differentiable.  Solving 
eq.  (3)  for  the  quadratic  error  criterion 


e(x-£.)  =  Tr{(x-^.)(x-^.)T] 


one  obtains 


de 

9Xi 


■2(x  -  £.) 


which  implies  that 


I  (x-£.)Px(x)dx  =  0 


i  =  1,2, ..  ,,M 


D. 

l 


Rearrangement  reveals 

J  2£P  (x)d2£ 

*i  = 


D. 


or 


!  P  (x)d* 
JD.  X 


^  =  E{x|xeDj 


i  =  1,2 . M 


(4) 


(5) 


(6) 


(7) 


(8) 


This  is  an  expression  for  the  best  nonlinear  mean  square  estimate  of  x, 
given  that  x  lies  within  region  D., 

Now  assume  that  x  is  distributed  according  to  a  Gaussian  probability 
density  function 


T  - 1 

p  (x)  =  K  exp[-  |x  C  x} 

X  ™  ™  ~  ”r 


(9) 


where  C  is  the  covariance  matrix  of  x  and  the  mean  is  assumed  to  be 
— x  — 

zero.  Also  let 


D.  =  {x.  |x.e [a.,  b.)] 
i  1  J  J  J 


j  =  1 , 2,  . . . ,  N 


(10) 


Then 
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(11) 


p— ■  t  - 1 

I  xK  exp{-jx  x}dx 


rT  ■  1 

K  exp{  -|x  x]  dx 


Curry  [2]  has  solved  this  equation  for  finely  quantized  values  of  x^.i.  e. 


b.  -  a.  <c. 
J  J  J 


j  =  1,2,.. .,N 


where  a.  is  the  standard  deviation  of  the  j  component  of  x.  His  approach 
J  “ 

is  to  approximate  the  Gaussian  density  by  the  first  three  terms  of  its  Taylor 

series  expansion  about  the  midpoint  of  the  region  D^.  The  integration  can 

then  be  performed,  with  the  result  that 

i  b+£. 

E{x|xeD.}  =  (I_-AC  )—  (13) 

1  ""“X  £ 


where 


(b.-a.) 

A  =  — J — ^ —  6 
-  12  kj 


k,j  =  1,2,. ..,N 


An  exact  solution  can  be  obtained  when  the  components  of  x  are  un¬ 
correlated.  In  this  case  the  covariance  matrix  can  be  expressed  as 


C  =  {of  6. .} 

-x  J  kj 


k,  j  =  1 , 2, . . . ,  N  (15) 


and  much  computation  reveals  that 


.  ..IT 

-i  V  TT 


al(e"bl/2al  -  e'al/2CTl) 
h~i  ~ 

erf  -erf  7^ 


2,„  2 


2  2 


a  (e"bN/2aN  .  e_aN/2aN) 
N _ 

f  bN  f 
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Gaussian  variables  which  had  been  decorrelated  by  means  of  a  Karhunen-Loeve 
transformation  and  then  quantized  could  be  restored  according  to  a  minimum 
mean  square  error  criterion  by  utilizing  this  last  equation. 

An  exact  analytical  solution  to  eq.  (11)  also  exists  when  an  estimate 
of  a  single  vector  component,  xN>  is  desired  based  upon  two  types  of  inform¬ 
ation  --  (a)  the  other  components,  x^x^...^  which  are  known 
completely  (quantized  with  an  infinite  number  of  bits);  (b)  the  quantizer  output 

which  nonlinear ly  specifies  the  interval  containing  x  .  To  derive  this,  consider 

N 


Nov/  denote  the  elements  of  (C  )_1  by 
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Then  performing  the  one -dimensional  integrations  in  eq.  (19)  yields 


■  exp  f- 


1  /  N"X  \2  I  bN 

ir~  rNNXN  +  y'  ViN  I 


NN  N 


Sj  is  Quantized  to  an  infinite  number  of  bits,  then  yr'1  =  a  =  b 
N  yi  N  N’ 

as  expected.  If  is  quantized  to  zero  bits,  its  interval  is  the  real  line 

(■aj^  =  =  “),  and  then  its  estimate,  y^t  is 

N  1 

yi  =  -7 —  y.  yjN  ( 

NN  “  J 
J  =  1 

This  result  is  identical  to  that  obtained  by  Pratt  [l]  in  estimating  an 
unknown  spectral  value  based  on  known  spectral  components.  However 
eq.  (20)  is  a  more  general  result  in  that  it  can  be  utilized  to  estimate 


components  that  have  been  quantized  to  any  number  of  bits  by  an  arbitrary 
quantization  scheme. 

Transform  Domain  Spectrum  Interpolation  The  above  solution  is 
applicable  to  the  mean  square  restc  ation  of  zonal  coded  transform  samples. 

In  this  case,  the  transform  samples  have  a  Gaussian  distribution,  since  each 
is  the  sum  of  a  large  number  of  random  variables  so  that  the  central  limit 
theorem  can  be  invoked.  These  transform  samples  are  typically  quantized 
according  to  a  bit  assignment  such  as  the  one  shown  in  Figure  1.  For  such 
a  quantizing  scheme,  only  eq.  (16)  can  be  utilized  directly  for  restoration; 
however  this  equation  ignores  the  known  correlation  existing  between  the 
samples.  Curry's  method  of  eq.  (13)  is  unable  to  restore  samples  quantized 
to  fewer  than  two  bits.  However,  for  greater  bit  assignments,  it  has  the 
advantage  of  providing  a  simultaneous  solution  utilizing  all  the  available 
information.  The  technique  developed  in  eqs.  (17)  to  (21)  avoids  the  above 
difficulties,  but  requires  a  recursive  solution  which  may  be  only  asymptotically 
optimal  (further  analysis  is  expected  to  establish  this).  Therefore  the  best 
restoration,  on  the  basis  of  optimality  and  ease  of  implementation,  is  obtained 
from  a  combination  of  the  solutions  presented  above  and  must  be  adapted  to 
the  particular  quantizer  used.  This  technique  will  soon  be  applied  to  zonal 
transform  coded  images.  It  is  anticipated  that  the  resultant  image  will  have 
a  lower  mean  square  error  and  improved  subjective  quality. 

References 

1.  W.  K.  Pratt,  "Transform  Image  Coding  Spectrum  Extrapolation, " 

Proc.  Seventh  Hawaii  International  Conference  on  System  Sciences, 
January,  1974,  pp.  7-9. 

2.  R.  E.  Curry,  Estimation  and  Control  with  Quantized  Measurements, 

The  M.I.T.  Press,  Cambridge,  Massachusetts  1970. 

3.6  Variable  Rate  Image  Coding  for  Sources  with  Unknown  Probabilities 
Lee  D.  Davisson 

The  average  distortion  of  image  encoding  at  a  fixed  rate  subject  to  a 


fidelity  criterion  depends  upon  an  actual  statistical  source  index  8,  in  effect, 
the  actual  stationary  ergodic  source  model  for  the  image  to  be  encoded.  Thus 
distortion  is  a  random  variable  over  the  ensemble  with  distribution  given  by 
the  distribution  of  the  parameter  8,  i.e.  the  class  of  possible  images.  In 
many  applications  it  may  be  more  desirable  to  allow  the  coding  rate  to  depend 
on  8  while  holding  the  average  distortion  fixed  over  the  ensemble.  This  is 
the  case,  for  example,  when  the  image  is  to  be  stored  or  a  variable  trans¬ 
mission  rate  exists  due  to  the  multiplexing  of  many  messages,  e.  g.  the 
ARPANET. 

A  coding  theorem  has  been  established  for  the  special  case  of  a  finite 

number  of  subsources,  8  =  1,  2, . . . ,  K,  in  the  ensemble.  The  theorem  holds 

for  noncountable  ensembles  as  well,  but  the  proof  is  considerably  more 

involved.  In  addition,  it  is  assumed  that  there  is  a  maximum  distortion 

value,  p  . 

M 

For  each  value  of  k,  generate  a  set  of  codewords  according  to  the 

usual  coding  theorem  for  stationary,  ergodic  sources.  If  R  (D)  is  the  rate 

th  ^ 

distortion  function  in  bits  of  the  k  subsource,  and  D  is  the  desired  value 

of  average  distortion,  each  code  will  contain  L,  codewords  where 

k 

log  Lk  =  (N(Rk(D-e)  +  e  )  (1) 

and  e  is  an  arbitrary  positive  constant  with  the  blocksize  chosen  large 
enough  so  that  the  average  distortion  is  D-e  for  all  8,  and  so  that  the 
probability  that  there  is  no  codeword  with  distortion  less  than  D-e/2  is 

less  than  e/2  p  . 

M 

The  coded  representation  of  each  of  the  L  codewords  generated  for 

k 

each  k  consists  of  two  parts.  The  first  part  is  the  fixed  length  binary 
number  equal  to  k-1,  k  =  1,  2, .  . .  ,  K  using  at  mo3t  log  K+l  bits  to  identify 
the  codeword.  The  second  part  is  the  location  of  the  codeword  in  a  list  for 
each  k  of  length  at  most  log  L^+l  bits.  Thus  the  rate  of  any  codeword  for  a 
given  0  is  at  most 
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Obviously  by  choosing  N  large  enough  and  c  small  enough,  the  rate  can  be 
made  arbitrarily  close  to  RQ(D)  for  0=  1,2,...,K. 

The  achievement  of  D  and  R0(D)  for  the  combined  supercode  depends 
upon  the  actual  choice  of  a  codeword  out  of  the 


total  codewords.  The  coding  rule  is  as  follows:  Upon  observing  an  output 
block  of  length  N,  among  the  codewords  of  distortion  less  than  D-e/2,  find 
the  one  of  minimum  rate,  if  any.  If  there  is  no  codeword  with  distortion 
less  than  D-e/2,  make  a  random  choice.  The  average  rate  for  0  -  k  then  is 


(3) 


which  is  arbitrarily  close  to  R^(D)  for  small  enough  e  and  large  enough  N. 

The  average  distortion  for  0  =  k  is 

Dk  SD-e/2  +  pM  Probjno  codeword  of  distortion  less  than  D  in  the  kth  code] 

DkSD-e/2+pMe/2pM  =  D  (4) 

The  result  described  is  largely  an  existence  theorem  and  therefore 
doer  not  prescribe  a  specific  method  of  synthesizing  data  compression 
systems.  It  does,  however,  provide  figures  of  merit  and  optimal  perform¬ 
ance  bounds  that  can  serve  as  an  absolute  yardstick  for  comparison  with 
real  systems.  Furthermore,  the  results  provide  theoretical  justification 
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for  overall  design  philosophies  that  have  proved  useful  in  practice.  In 
practice,  one  could,  for  example  have  a  coding  scheme  for  each  of  several 
classes  of  pictures,  e.  g.  classed  according  to  "busyness.  "  The  appropriate 
coder  would  be  switched  in  for  each  encoded  image  with  the  coder  identity 
sent  as  a  prefix  to  the  encoded  picture  information  as  suggested  above. 

The  approximate  minimum  number  of  bits  would  then  be  used  for  each 
image  depending  upon  the  average  allowed  distortion. 


Image  Restoration  and  Enhancement 

Image  restoration  techniques  seek  to  reconstruct  or  recreate  an 
image  to  the  form  it  would  have  had  if  it  had  not  been  degraded  by  some 
physical  imaging  system.  Image  enhancement  techniques  have  two  major 
purposes:  improvement  in  the  visual  quality  of  a  picture  to  a  human  viewer; 
and  manipulation  of  a  picture  for  more  efficient  processing  and  data  extrac¬ 
tion  by  a  machine.  Both  techniques  are  subjects  of  continuing  study;  re¬ 
sults  of  this  effort  during  the  pa3t  six  months  are  summarized  below. 

The  first  report  deals  v  ith  methods  of  restoration  utilizing  the 
matrix  pseudoinverse  of  a  blur  matrix  which  models  a  space  variant 
point  spread  function.  Computation  of  the  pseudoinverse  is  obtained  indir¬ 
ectly  by  a  singular  value  decomposition  of  the  blur  matrix.  The  following 
report  considers  another  approach  to  pseudoinverse  image  restoration.  Com¬ 
putational  techniques  are  developed  for  pseudoinverse  restoration  by  process¬ 
ing  in  the  transform  domain.  « 

An  anal  sis  is  presented  in  the  next  report  of  the  effects  of  discrete 
modelling  of  the  superposition  integral  for  image  restoration.  The  effect 
of  modelling  errors  and  the  ill- conditioning  of  the  blur  matrix  are  quan¬ 
tified  for  typical  image  blur  models. 

The  use  of  spline  functions  in  imaging  models  is  explored  next.  It 
is  shown  that  the  B- splines  offer  computational  promise  for  image  restor¬ 
ation  with  smoothness  constraints. 

Astigmatism  and  curvature  of  field  image  degradations  are  character¬ 
ized  by  spatially  variant  imaging  models.  It  is  shown  that  the  imaging  model 
can  be  decomposed  into  a  cascade  of  a  geometric  distortion,  a  spatially 
invariant  linear  system,  and  another  geometric  distortion.  Methods  of  in¬ 
verting  this  cascaded  model  are  developed. 

In  the  next  report,  the  concept  of  transform  domain  Wiener  filtering 
for  image  restoration  is  extended  to  include  lower  triangular  transformations. 
For  a  Markov  process  image  model,  it  is  found  that  an  extremely  efficient 
computational  algorithm  can  be  obtained. 

The  two  following  reports  are  concerned  with  two  aspects  of  color 
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image  restoration.  Once  report  considers  techniques  for  estimating  the 
tristimulus  values  of  a  color  from  spectral  observations  of  the  color.  The 
other  report  presents  methods  for  compensating  for  film  nonlinearities  in 
a  photograph  of  a  color  television  display. 

The  final  report  presents  a  quantitative  development  of  pseudo¬ 
color  image  enhancement  techniques  to  map  a  grey  scale  image  into  a 
color  display.  Mappings  are  found  that  are  psychophysically  relevant,  yet 
computationally  efficient. 


4.1  Space  Variant  Point  Spread  Function  Pseudoinversion 
Monty  Adler  and  Harry  C.  Andrews 

Many  complex  forms  of  image  degradation  cannot  be  modelled  by 
a  spatially  invariant  point  spread  function;  consequently,  Fourier  tech¬ 
niques  are  not  applicable  to  the  restoration  process.  For  these  spatially 
variant  point  spread  function  systems,  restorations  can  be  achieved  by  a 
matrix  formulation  in  which  point  spread  function  matrices  are  inverted. 
For  singular  blur  matrices,  the  inverses  must  be  replaced  by  pseudoin¬ 
verses  to  achieve  the  least  squares  approximation  to  the  original.  For 
computational  simplification  an  assumption  of  a  separable  space  variant 
point  spread  function  allows  the  following  analysis. 

Let  F  be  a  matrix  representing  a  perfect  image  which  is  acted  upon 
by  a  separable  blur  function  to  produce  a  blurred  image  (j  as  modelled  by 

G  =  A  F  B 

where  13  blurs  the  rows  of  F  and  A  blurs  the  columns  of  F.  Usually  A  and 
B.  are  singular  or  almost  singular.  One  would  like  to  find  an  approximation 

A  A 

to  F,  F  such  that  ||F  -  F  ||is  minimized.  The  analysis  and  results  which 
follow  do  not  depend  on  any  of  the  matrices  being  square  but  for  simplicity 
it  will  be  assumed  that  all  matrices  are  N  by  N.  If  A  and  B  are  nonsingular 
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A  -1  -1 

and  square  the  solution  is,  of  course,  F  =  A  G  B  =  F. 

Using  appropriate  A  and  B  matrices  one  can  model  space  variant 
blur,  linear  motion  blur,  and  in  general  any  separable  blur  function. 

It  is  known  that  A  and  B  can  be  expressed  as  a  sum  of  matrices  of 
rank  one  as  follows 


N 

A  =  E  ^  Ua  VaT 
~  i  =  1  1  ^  ^ 


for  \  *  Vi  2  ° 


B  :  ]T  rf  UbVbT 
"  i  =  l  1  -1"1 


,  b  b 

for  r\.  2  T). ,  a  0 

i  l+l 


where  TL  ,  Tib  are  scalars,  and  uf  ,  Ub,  V.a,  V.b  are  column  vectors  of 

11  a  h — 1  — 1  1  W"""1 

length  N,  and  the  sets  [uf]  ,  {U.  }  ,  £va}  ,  {VT)are  orthonormal.  In  an  ideal 

— i  — l  “i  — i 

computational  environment,  if  =  0  for  i  >  R  where  R  is  the  rank  of  A; 

i  a  a  — 

similarly  for  B. 

The  ideal  pseudoinverse  of  A  can  be  expressed  as 


N  + 

a+  =  £  n.a  vf  ua 

iTi  ‘ 


where 


a+  1,  . ,  a  . 

n.  =  /  a  if  ri.  t  0 

1  Tl.  1 

1 

=0  if  T| a  =  0 

i  i 


A  model  is  needed  for  the  computed  SYD.  Using  a  computer  the  calculated 
SVD  will  be  modeled  by 


N 


A  =  E  Xa  Ua  for  £  2  £  2  0. 

.  ,  i  —i  — i  i  i+I 

i  =  l 


£1  A 

The  X  will  not  be  zero  for  i  >R^due  to  computational  inaccuracy;  and  Uf  , 
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vf  will  in  general  be  identical  to  the  ideal  case.  Note  that  the  rank  of  A , 

2i 

R>a  cannot  be  accurately  determined  since  the  X.  do  not  go  to  zero.  Defining 


.  + 
—  K 


K 


y  i/xa  va  ua 

i  =  i  1  -1  -1 


aT 


then  as  soon  as  K  goes  beyond  R  ,  the  algorithm  divides  by  1/X.  for  an 

■  a  1 
&  T 

inaccurate  X  and  A  blows  up.  It  is  desirable  to  "  stop  K  at  R  "  but  R 
i  is.  a  a 

is  not  known. 


The  problem  is  further  complicated  if  A  and  15  are  non-singular 

&  fo 

but  have  small  X^  ,  X^  which  cannot  he  computed  accurately.  Even  if  R^  and 

R,  were  known,  the  algorithm  would  not  want  to  divide  by  the  inaccurate 
a  b 

X.  ,  X.  values.  Unfortunately  the  typical  situation  is  that  A  and  B  are  either 

a  b 

singular  or  non- singular  with  small  X^  ,  X  ;  but,  of  course,  it  is  not  known 
which  is  the  case. 

For  simplicity  A  and  B  will  be  assumed  to  have  the  same  or  close 
to  the  same  rank.  The  methods  described  below  are  trivially  modifieable 
if  the  ranks  of  A  and  B  differ. 

Define 


and 


£jk  *  4  ^  4 


2jk  *  ASKS  *  -A^J-Ci-B 


Since  under  the  above  assumption  the  best  J  equals  the  best  K  it  is  possible  to 
drop  the  double  subscript. 


V  -K-^K 


G  =  A  F  B 
-K  - -K  - 


It  is  desirable  to  minimize  ||  F  -  F  ]]  over  all  values  of  K.  Three  methods 

- tv 
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are  proposed,  one  of  which  can  be  shov/n  to  find  the  best  K. 

A  A 

Method  1  -  Human  Intervention.  Successively  compute  F  ,  F  ... 

*  —  1  —  2 

and  display  to  a  human.  Although  the  best  F  will  not  be  apparent,  it 

A  I\.  A 

will  be  clear  when  F  is  blowing  up.  In  figure  1,  the  estimate  F  is 

R  51 

seen  to  blow  up,  and  this  fact  is  obvious  to  a  viewer.  For  a  display 

purposes  all  F^  values  have  been  rounded  to  the  nearest  integers,  all 

negative  values  have  been  set  to  zero,  and  all  values  greater  than  63 

have  been  set  to  63.  This  explains  the  black  and  white  blocks  in  the  figures, 

which  are  out  of  range  values. 

Method  2  -  Track  A  A^  and  JB^  B.  It  can  be  shown  that 


Pi  -  A  A^.  ||  2  =  N  -  K 


where  N  is  the  size  of  A.  Then  one  would  expect  E  (A)  =  ||l-A  A+  jj2  -(N-K) 

to  deviate  from  zero  when  A^  blows  up.  Trials  on  a  computer  show  that 

remains  at  a  low  constant  value  until  K  approaches  R  at  which  point 

increases  by  a  factor  of  10  around  R  .  It  thus  appears  that  tracking 

+  +  a  6 

—  —  or  — k  —  could  be  a  method  of  find  the  best  K.  Table  1  contains  the 

results  of  two  a  typical  computer  trial  where 

E(A)  =  1010  (l|l_- A^_A  ||  2  -  (N-K)  ) 

E(B)  =  1010  (III-  B+kB  ||2  -  (N-  K)  ) 

E(F)  =  ||_F  -  Fr  || 

E(G)  =  ||  G  -  Gk  || 

A 

Method  3.  Track  j|  — K  jj  It  can  be  shown  that  min  ||_F  -  F^  || 

occurs  at  the  same  value  of  K  as  min  ||  G  -  G  ||  and  therefore  the  best  K 

K 

can  be  found  simply  by  tracking  min  ||  G  -  G  l| .  Figures  1  to  3  show  the 

K  ^ 

results  of  three  simulation  cases.  In  each  case  the  best  reconstruction 
occurs  at  the  K  for  which  ||G  -  ||  is  a  minimum. 

X\. 
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(e)  (best  restoration)  (f)  F-.  (note  space  variant 

singularity) 

Figure  4.1-1.  Restoration  of  space  variant  point  sources,  N=60 

(Gaussian  blur  best  at  center)  (Computer  Generated  Image). 


I 


(e)  F^(best  restoration) 


(f)  F^g(note  space  variant 
singularity) 


Figure  4.  1-2,  Restoration  of  spa  variant  point  sources,  N=60 
(Gaussian  blur  best  at  center). 
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(e)  (f)  F^Q2(minimum  error) 


Figure  4.  1-3.  Restoration  from  space  variant  point  spread  function,  N-l  12 
(blur  best  at  center) 


Table  I  Paeudoinverse  and  Reconstruction  Computation  Error  - 


I 

E(F) 

E(G) 

E(A) 

E(B) 

0 

156.9 

6 

166.  3 

25.  6 

2.5 

2.4 

12 

158.9 

12.  8 

2.2 

2.  2 

18 

147.  2 

1.  9 

2.  0 

2.  0 

24 

141.6 

0.  3 

1.  7 

1.  7 

30 

134.  8 

2.4E-2 

1.  3 

1.  3 

36 

125.9 

1.1E-3 

1.1 

1.1 

42 

114.0 

1.6E-5 

1.2 

0.13 

43 

109.7 

7.  6E-6 

-0.1 

2.  5 

44 

108.4 

2.  2E-6 

-4.5 

2.  9 

45 

106.  5 

1.  6E-6 

6.  3 

-3.8 

46 

104.  3 

3.1E-7 

8.8 

-16. 

47 

101.1 

2.  589E-7 

58. 

-33. 

48* 

99.  2 

2.  583E-7 

-100. 

94. 

49 

177.  8 

1.  9E-6 

-67. 

510. 

50 

382.  9 

2.  9E-6 

-420. 

-250. 

51 

2.  2Et4 

2.  3E-  5 

-8800. 

7200. 

*  minimum 

Experimental  Results  Figures  1,  2,  and  3  show  an  image,  its  blur 
and  the  reconstruction  for  different  values  of  K.  Note  that  once  the  minimum 
is  found,  the  reconstructed  picture  quickly  blows  up.  Table  1  shows  ||  G  -  G  | 
and  ||  F  -  F^  ||  for  different  K.  Although  the  minimum  occurs  at  the  same 
point,  ||  G  -  Gr  ||  remains  small  after  the  minimum  by  ||  F  -  F  ||  blows  up. 
Naturally,  in  the  real  world  one  would  not  be  able  to  track  ||  F  -  F  ||  but 
II  G  -  G  ||  is  available. 


4.2  Pseudoinverse  Image  Restoration  by  Transform  Domain  Processing 
William  K.  Pratt 


Linear  operations  on  data  can  often  be  performed  more  efficiently 
by  indirect  techniques  which  involve  projection  of  the  data  to  another  vector 
space  through  a  unitary  transformation  of  the  data.  This  concept  has  been 
applied  with  success  for  discrete  Wiener  filtering  of  images.  The  following 
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outlines  the  extension  of  the  concept  to  pseudoinversion  for  image  restoration 
Imaging  Model  Let  F(x,y)  denote  an  ideal,  continuous,  infinite  ex¬ 
tent  image  field,  and  let  F  (x,  y)  represent  an  observed  image  field  which  is 
also  continuous  and  of  infinite  extent.  For  a  large  class  of  imaging  systems 
the  observed  image  is  related  to  the  ideal  image  by  the  convolution  integral 
~  00  00 

F:x,y)  =  J  J*  F(a,  B)g(x-a,  y- (3  )dadB 


where  g(x-a;y-g)  represents  the  impulse  response  of  the  imaging  system. 

In  the  discrete  model  of  the  imaging  system,  the  observed  image  is  repre¬ 
sented  by  physical  samples  spaced  evenly  over  a  unit  grid,  and  the  continu¬ 
ous  integration  is  approximated  by  a  quadrature  formula  resulting  in 

L+rrij-1  L+m^-l 

F(mi'm2)  =  £  £  F(n  ,n  )H(m  -n  +L,m.-n,+L)  (1) 

n^rrij  n2=m2  11  22 

fol  and  lsn.=N  where  the  array  H,  assumed  to  be  aero  outside  its 

rang,  of  indices,  represents  the  sampled  impulse  response  and  incorporates 
all  quadrature  factors.  The  impulse  response  is  also  truncated  to  an  LxL 

array.  In  order  to  prevent  serious  modelling  errors  at  the  image  boundary, 
it  is  necessary  that 


N  &  M  +  L  -  1 


It  is  computationally  convenient  to  represent  the  data  arrays  F  and  F 
as  column  vectors,  f  andj,  respectively  by  column  scanning.  Then  the  sam 
pied  superposition  operation  can  be  described  by 


where  B  is  an  M  X  N2  matrix  which  can  be  partitioned  as 
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0 


0 


0 


B  = 


-1,1  -1,2  —1,  L 

0  -2,2  -2,L  B2,  Lrf-1  °  •**  0 


0  0 


—  M,N 


(3) 


The  general  term  of  13  is  then  given  by 


Bm1.n1(m2'n2>=H,ml-VI'-m2-VL) 


(4) 


for 


1  s  ^  M  1  ^  s  M 

^  nj  s  L+m^  -1  s  n2  ^  L+m^  _1 

Pseudoinverse  The  imaging  model  of  eq.  (2)  can  be  "inverted" 

+  A 

by  a  pseudoinverse  operator  13  in  the  sense  that  an  estimate  {_  of  the  ideal 
image  vector_f  can  be  computed  by 

£  =  B+  l_  =  B+  B  (5) 

2 

If  13  is  of  rank  M  ,  then  the  pseudoinverse  is  equal  to 

+  T  T  -1 

B  =  B  (B  B  )  (6) 

2  2 

It  should  be  noted  that  since  M  <N  ,  that  is,  there  are  fewer  observations 
than  points  on  the  ideal  image  to  be  estimated,  the  estimate  will  not  be  exact 
even  in  the  absence  of  observational  error.  The  pseudoinverse,  however, 
does  provide  a  minimum  mean  square  error,  minimum  norm  estimate. 

Transform  Domain  Pseudoinverse  Figure  1  illustrates  the  compu¬ 
tational  steps  involved  in  direct  pseudoinversion  of  an  image  vector,  and  in 


B+ 

o)  IMAGING  MODEL 


A 

B  + 


b)  DIRECT  PSEUDOINVERSE  PROCESSING 


c)  TRANSFORM  DOMAIN 

PSEUDOINVERSE  PROCESSING 


Figure  4.2.1.  Direct  and  tran.form  domain  pseudoinver.e  proceaeing. 
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the  transformation  processing  approach.  In  transform  processing  a  uni¬ 
tary  transformation  i'j  performed  on  the  observed  vector  £  prior  to  multi¬ 
plication  by  the  transform  domain  pseudoinverse  matrix  An  inverse 

A 

transform  reconstructs  f.  From  figure  1,  since 


f  =  r 


and 


then  clearly 


i  •  (AmZ)'1  <An2>  i 


1  =  (A  2)B+  (A,) 


-1 


M 


N 


It  is  easy  to  show  that 


6+  =  6*T  B*T) 


where 


£  =  (A  B  (A  ,) 

~UC - 


-1 


(7a) 

(7b) 

(8) 


(9) 

(10) 


is  the  transform  domain  representation  of  the  blur  matrix  of  eq.  (3). 

Computational  efficiencies  in  transform  domain  pseudoinversion 
result  from  the  sparseness  and  structure  of  the  blur  matrix  B.  As  an 
example,  consider  the  Fourier  transform  represfentation  of  I3+.  In  this 
case  the  transformation  matrix  is  of  the  form 


where 


®  — K 


(ID 


=  _L  w(x-l)(y-l) 
K 


r(N) 


for  x,  y  =  1,2,...,  K.  Now,  let  Hj,  denote  the  extended  impulse  response 
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PM 


obtained  by  imbedding  the  L  X  L  impulse  response  matrix  in  the  upper  left 
corner  of  an  N  X  N  matrix  of  -.eros.  The  two  dimensional  Fourier  transform 
of  the  extended  impulse  response  matrix  is  obtained  from 


V<N> 

-E 


A  H(N) 

-N  — E 


(12) 


These  transform  components  are  then  column  scanned  and  inserted  as  the 
diagonal  elements  of  the  N2  x  N2  matrix 

v(N)  ..  rl/ 

-  =  dia*  CE(M).ifE(2,l) . Ve(N,N)}  (13) 

Then,  it  can  be  shown,  after  considerable  manipulation,  that  the  Fourier 
transform  blur  matrix  is 


B  =  tB®PB]v(N> 


(14) 


where 


J?B(U*  v) 


,  w  -(v-l)(L-l) 
1  A  "  WN 


Vm  i  -  w  (u_1)w  ■(v_1) 
M  N 


(15) 


Thus,  the matrix  operator  consists  of  a  scalar  weighting  matrix,  V(N), 
and  a  matrix  CPg  ®  Pg]  that  performs  the  dimensionality  conversion  (an 
interpolation  operation)  between  an  N2  element  input  vector  and  an  M2  eleme  nt 
output  vector.  The  dimensionality  matrix  is  extremely  sparse,  and  therefore, 
savings  can  be  obtained  in  the  computation  of  £  and  subsequently,  _£+.  As  an 
example  figure  2  contains  displays  of  the  blur  matrix  B,  the  pseudoinverse 
matrix  B  ,  and  their  transform  domain  representations,  B_  and  -9+  for  a  9  x  9 
Gaussian  shaped  impulse  response,  and  M=8  and  N=l6.  The  relative  sparse¬ 
ness  of  6  and  8  are  apparent  from  the  figures. 
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(a)  spatial  domain, 


(b)  Fourier  domain,  B 


\ 


(c)  spatial  domain, 


(d)  Fourier  domain,  B 


1 


Figure  4.  2-2.  Examples  of  blur  matrix  and  its  pseudoinverse 
in  spatial  and  Fourier  transform  domain. 


4.3  Modelling  Superposition  Integrals  for  Image  Restoration 
Nelson  D.  A.  Mascarenhas 

The  process  of  image  blurring  and  addition  of  noise  in  an  incoherent 
optical  system  can  often  be  described  by  the  equation 

y(a*  P)  =  J_b  J,b  x(?,ri)h(af?;ftT))d?dTi+Ti(af0  -®<a,  p<® 

U  U 


where  x(a,  P)  denotes  an  ideal  image,  y(a,  P)  is  the  observed  image,  p  (a,  P) 
is  an  additive  noise  field,  and  h(a,  ?;P,Ti)  is  the  image  system  spread  function. 
When  a  restoration  problem  is  to  be  solved  with  a  digital  computer  to 
estimate  x(a,  P),  a  discretization  has  to  be  performed.  By  using  a  lexico¬ 
graphic  ordering,  it  is  possible  to  red-ce  the  resulting  two  dimensional 
data  arrays  into  vector  format.  The  following  equation  describes  the  dis¬ 
crete  model 


_y  =  B  x  +  n 

where 

2 

_y  =  (M  x  1)  vector  of  observed  values 

„  2  2 

13  =  (M  X  N  )  blur  matrix 

2 

x  -  (N  XI)  vector  of  original  pixel  values 
2 

n  =  (M  XI)  vector  of  noise  components. 

In  general  the  entries  of  the  blur  matrix  depend  on  both  the  kernel  of  the 
integral  equation  and  the  weights  of  quadrature  integration.  In  the  simula¬ 
tion  experiments  described  in  this  section  these  weights  have  been  assumed 
to  have  unity  value. 

Figure  1  describes  the  data  arrays  involved  when  an  overdetermined 
model  for  restoration  is  used.  This  leads  to  the  following  structure  for  the 
blur  matrix  B 
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ORIGINAL  PICTURE 


Figure  4 


> 


3-1.  Data  arrays  in  the  overdetermined  model 
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First  B  is  partitioned  into  submatrices  .  of  size  (M  X  N).  Then  each 
submatrix  has  a  similar  structure,  being  composed  by  a  nonzero  diagonal 
band  of  elements. 

Two  expressions  for  the  blur  have  been  used.  The  first  simulates 
the  effect  of  atmospheric  turbulence  over  a  long  exposure.  The  spread 
function  is  given  by 

h(al'?m:Pj'V  =  6XP  ‘“NT  +  ~ NT- 


where  the  coefficients  b„  and  b„  control  the  amount  of  blur  imposed  on  the 

V  ri 

vertical  and  horizontal  directions,  respectively.  The  second  blur  function, 
also  space  invariant  and  in  separable  form,  simulates  the  effect  of  a  diffrac- 
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For  white  noise,  the  bent  linear  unbiased  estimator  of  the  original  picture 
is  obtained  as 

A  ._T_.-1._T 
x  =  (B  B)  B 

and  the  covariance  matrix  of  this  estimator  is 

2  T  -1 
Va  =  a  (B  B) 

2 

where  a  is  the  variance  of  the  noise.  The  amount  of  uncertainty  on  the 

T 

estimators  will  depend  on  the  degree  of  singularity  of  (B  B).  A 
possible  measure  of  this  is  given  by  the  condition  number  of  the  matrix _B, 
which  can  be  expressed  by  the  ratio  of  the  largest  to  the  smallest  singular 
value  of  the  matrix  15  [l]. 

Figures  2  and  3  illustrate  curves  of  condition  number  versus  blur 
coefficient,  for  a  given  number  of  original  pixel  values  (N=8)  [2].  In  the 
curves  the  number  of  observed  values  M  is  varied  while  maintaining  the 
structure  of  the  blur  matrix  described  previously.  The  existence  of  a 
maximum  of  the  condition  number  curves  can  be  explained  in  terms  of  the 
truncation  of  the  point  spread  function  displayed  in  Figure  4.  In  fact,  for 
increasing  M,  the  number  of  points  where  this  function  can  be  nonzero  is 
increased,  and  the  effect  of  the  truncation  starts  only  for  higher  blur  coef¬ 
ficients.  Consequently,  the  curves  for  different  values  of  M  have  essentially 
a  common  ascending  branch  and  the  descending  part  starts  at  varying  parts 
for  different  values  of  blur  coefficients.  If  there  were  no  truncation,  the 
curve  would  approach  infinity  very  fast,  the  asymptotic  value  being  obtained 
for  the  smoothest  possible  kernel,  with  constant  value  unity,  implying  a 
blur  matrix  with  rank  one.  With  the  truncation,  the  curves  show  a  descending 
branch  that  begins  at  the  point  where  the  increasingly  wider  kernel  starts  to 
be  cut  down  substantially.  For  greater  amounts  of  blur,  the  curves  tend  to 
a  finite  value  that  depends  on  M. 

The  curves  of  figs.  2  and  3  can  be  used  as  a  guide  for  the  choice  of 
the  number  of  sampling  points,  once  the  number  of  quadrature  nodes  is  fixed. 


! 


Figure  4.  3-2.  Condition  number  curves  for  different  number  of  sampled 
values  -  Gaussian  blur* 
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jj* 


CONDITION  NUMBER 


spread  function. 


For  a  very  small  amount  of  blur  all  curves  coincide  so  that  the  designer 
may  choose  M  =  N  with  almost  no  error.  In  this  case  blur  plays  no  role, 
only  noise  affects  the  restoration.  With  increasing  blur,  different  numbers 
of  sampling  points  give  different  values  of  condition  number.  If  a  curve  on 
an  ascending  branch  is  chosen,  truncation  has  no  effect  on  the  kernel,  but  a 
high  condition  number  imposes  high  variances  on  the  estimators.  If  a  curve 
on  a  descending  branch  is  selected,  lower  variances  of  the  estimators  are 
obtained,  at  the  price  of  error  on  the  estimation  of  the  continuous  function 
due  to  the  truncation  error  in  the  discrete  model.  Therefore,  a  trade-off 
between  the  variance  of  the  estimators  and  the  modeling  error  can  be  char¬ 
acterized. 

Although  these  conclusions  are  drawn  based  on  the  particular  model 
discussed  in  this  section,  they  are  more  general.  Since  the  inverse  of  the 
integral  operator  that  describes  the  blur  is  unbounded,  therefore,  the  closer 
the  discrete  model  follows  the  continuous  one,  the  more  ill  conditioned  the 
former  model  tends  to  be.  A  move  in  the  opposite  direction  reduces  singu¬ 
larity  but  imposes  modeling  errors.  This  inevitable  dilemma  can  only  be 
broken  with  correct  a  priori  knowledge  about  the  solution. 
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4.4  Spline  Function  Restorations 

Steve  Hou  and  Harry  C,  Andrews 

A  variety  of  models  for  linear  imaging  systems  have  been  postulated 
for  processing  by  a  general  purpose  digital  computer.  Three  of  these  models 
are: 
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Continuous -continuous.  Let 


g(*»  y)  =  Ti)h{x,  y,  C,  n)dCd  n 

where  the  image,  g,  object,  f,  and  point  spread  functions,  h,  are  all  des¬ 
cribed  in  continuous  notation. 

Pis  Crete -discrete  .  Let 


£  =  Hi 

where  the  image  and  object  have  been  scanned  or  stacked  into  x  1 

2  2 

column  vectors  and  H  is  an  N  x  N  matrix. 

Continuous-discrete.  Let 

9.  =  H  H(  C#T))f(  C,h)dC  dp 

where  the  image  is  a  matrix  of  entries  g..  and  the  object  is  continuous  and 
passes  through  a  matrix  of  continuous  point  sprerd  functions.  Jhus 

gij  =  JJ^j^WMdCdri 


Figure  1  indicates  the  geometry  associated  with  this  model.  Possibly 
this  model  is  the  most  realistic  of  the  above  alternatives  in  the  sense  that 
the  object  is  continuous  before  being  sampled  and  the  image  (once  in  the 
computer)  is  inherently  discrete. 

In  attempting  a  restoration  of  the  discrete  image  G  to  a  better 
estimate  f(C,  n)  of  the  continuous  object  f(C,  n)  it  is  clear  that  the  degrees 
of  freedom  in  the  object  are  infinite  while  those  in  the  image  are  finite. 
Thus  it  is  necessary  to  provide  a  translation  of  the  approximation  of  a 
continuous  function  to  the  estimation  of  a  finite  representation  of  a  func¬ 
tional  form  of  that  continuous  object.  The  mechanism  of  spline  functions 
are  suggested  for  this  purpose.  If  f(C,n)  denotes  the  restored  object,  then 


let 


«•  « 


p  =  JJflrc,Ti)f(C,n)dCdTi 


It  is  then  possible  to  formulate  the  restoration  expression  in  terms  of  the 
minimization  of  an  objective  function 

W(f)  =  tr{[G-P]T[G-P]  } 

In  this  case,  the  estimation  is  a  least  squares  restoration.  Additional 
constraints  on  the  objective  function  may  be  desired.  One  such  constraint 
might  be  that  the  restored  object  not  oscillate  wildlv;  in  particular  that  the 
second  derivative  of  that  object  be  minimum.  Then  using  the  method  of 
Lagrange  multipliers  the  objective  function  becomes 

W(f)  =  tr{fe-P]TCG-P]}+>J‘J,(V2f(Cti)  )2dCdri 

Now  expanding  the  restored  object  into  a  set  of  nonorthogonal  basis  func¬ 
tions  determined  by  bi-cubic  B  splines  on  a  uniform  two  dimensional  grid, 
one  obtains 


f(C,  -n) 


caeSa(C>¥"> 


where  is  the  cubic  B  spline.  Such  an  expansion  allows  a  priori  know- 

* 


ledge  of  the  second  derivative  of  f  (£,  n)  in  terms  of  the  coefficients  c 
When  the  objective  function  is  minimized,  the  matrix  equation 


ap- 


(A  +  lB)c  =  d 

is  obtained  where  A_  is  defined  by  the  blurred  spline  functions,  B  is  the 
differential  spline  functions  and  highly  banded,  £  is  the  unknown  coefficient 
vector  and  <3  is  defined  by  the  image  G. 

Techniques  for  solving  c,  the  equality  constrained  least  squares  prob¬ 
lem  described  above,  are  being  investigated  in  the  framework  of  spline  func¬ 
tions  and  image  arrays.  In  addition  inequality  constraints  on  the  c  vector 
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such  that  f(C,  n)  is  positive  and  bounded  are  being  considered.  The  approach 
looks  fruitful  as  the  B- spline  basis  functions  provide  a  unique  Met  of  non- 
orthogonal  basis  functions  which  are  inherently  positive  and  whose  coeffi¬ 
cients  provide  second  derivative  (smoothness)  criterion  control. 

4.  5  Space-Variant  Restoration  of  Astigmatism  and  Curvature  of  Field 
Alexander  A.  Sawchuk  and  M.  Javad  Peyrovian 

As  discussed  in  previous  work,  the  general  problem  of  image  restor¬ 
ation  with  a  space-variant  point- spread  function  (SVPSF)  deg  adation  is 
exceedingly  difficult  to  solve  due  to  the  high  system  dimer sionality  &,  2], 

To  get  a  practical  solution,  it  is  necessary  to  completely  exploit  the  de¬ 
grading  system  symmetry  and  structure  in  order  to  reduce  the  effective 
dimensionality.  The  approach  has  wide  applications  in  the  restoration  of 
motion  blur  &,  3]  and  certain  geometrical  optics  aberrations  \Z,  5,  61 .  This 
section  describes  an  extension  of  the  technique  to  astigmatism  and  curvature 
of  field. 

A  great  deal  of  insight  into  the  restoration  problem  is  obtained  by 
careful  examination  and  derivation  of  the  degrading  system  point- spread 
function  (PSF).  The  aberration  function  &,  2,  4l  of  geometrical  optics  for 
astigmatism  and  curvature  of  field  are  given  by 

VU1  =  A1  =  (2C+D)u^  ercos  (la) 

VU2  =  *2  =  D  "l  S  Sin  E0  (lb) 

where  and  (u^,  u^)  are  the  rectangular  image  and  object  coordinates 

respectively,  and  are  ray  intercepts  in  the  exit  pupil  of  the  optical 
system,  and  C  and  D  are  constant  coefficients  describing  the  degree  of  astig¬ 
matism  and  curvature  of  field,  respectively.  Using  previous  results  Q,2l 
an  equivalent  system  SVPSF  h(x,  u]  which  describes  the  degradation  can  be 
derived.  For  a  circular  exit  pupil  of  radius  R  this  SVPSF  is 


-6  8. 


h(xj,  x2>  u  ,u2=  0)  = 


0(20+0)^ 

0 


R2R2u^ 


2  4  ? 

(2C+D)  ”  R^ 


*  1 


elsewhere 


(2) 


with  the  input  object  impulse  function  at  (u  ,u2=0).  The  function  of  eq.  (2) 

is  shown  in  figure  la  for  the  impulses  at  various  locations  in  the  (u  u  1 

1  2 

plane.  The  regions  of  nonzero  response  are  defined  by  ellipses  which  in¬ 
crease  in  size  as  the  square  of  the  radial  distance.  At  the  same  time, 
the  amplitude  of  the  response  decreases  inversely  with  u* .  Because  of 
the  inherent  circular  symmetry,  the  amplitude  and  shape  of  the  response 
is  a  function  of  radial  distance  only  and  not  a  function  of  angle  8,  as  shown 
in  figures  la  and  lb.  Thus  one  reduction  of  system  complexity  is  obtained 
bY  rewriting  the  aberration  functions  of  eq.  (1)  in  polar  coordinate  form  as 


x  -u  =  (2C+D)u  e  cos  en 
r  r  r  r  8 


(3a) 


-1 


XQ-U0  =  *an  (Du^e^sin  Eq/I+^C+DJu^e^cos  Eg)  (3b) 


where  (xr*xg)  and  (ur,Ug)  are  image  and  object  polar  coordinate  variables. 
This  procedure  suggests  that  the  coordinate  transformation  restoration  (CTR) 
techniques  used  for  other  kinds  of  space- variant  degradations  are  also  valu¬ 
able  here.  Performing  this  transformation  to  the  system  described  by  eq.  (3) 
effectively  converts  the  problem  to  a  space-invariant  blur  in  8  which  changes 
slowly  with  u^(eq.  (3b)  )  and  a  purely  radial  two-dimensional  space-variant 
blur  (eq.  (3a)  ). 

For  the  case  of  pure  curvature  of  field  with  no  astigmatism,  C  =  0  in 
eq.  (2)  becomes  singular.  The  SVPSF  for  astigmatism  only  has  no  8  blurring 
and  collapses  to  a  radial  space-variant  blur  h^,  u^,  u2=0)  obtained  by  evalu¬ 
ating 


u. 

hc(xruru2 =0)  =  j  h(xrWu2  s  0)  dxj 


(4) 


and  letting  D  approach  0.  The  result  is 


hc(xr  ui’  u2=0)  = 


9 


(4C  R  U1  '  (X1'U1}  ^  u^ZCu*  R  ui+2Cuj  R 


9r2  4 
ZC  Uj 


0 


elsewhere 


(5) 


which  is  shown  in  figures  2a  and  lc.  For  astigmatism  only,  there  is  no 
blurring  in  the  A  direction  anu  one-dimensional  space- variant  blur  in  the 
radial  direction. 

The  reduced  dimensionality  of  the  astigmatism  degradation  following 
the  polar  coordinate  transformation  makes  it  possible  to  write  the  degrada¬ 
tion  operation  in  matrix  form  as 


Iix9»xr)  =  ii(xr*url2.(x0»ur)  (6) 

wh«*r«i  i(xQf  x  )  and^fx^u  )  are  one-dimensional  line  images  for  each  xn  and 
Hfx^.u^)  is  the  discrete  blur  matrix  obtained  applying  a  quadrature  formula  or 
other  approximation  to  the  continuous  space  description  of  eq.  (5).  The  degra¬ 
dation  due  to  astigmatism  and  the  restoration  process  is  shown  in  a  series  of 
figures  that  follow.  Figure  2b  shows  an  original  undegraded  object  scene, 
and  fig.  2c  shows  the  effect  on  the  original  of  space-variant  astigmatism 
described  by  eq.  (5).  Note  that  the  amount  of  blurring  increases  proportional 
to  the  square  of  the  distance  from  the  origin. 

One  method  of  image  restoration  for  the  space- variant  system  of  eq.  (6) 
is  to  use  a  discrete  pseudo-inversion  technique.  As  a  result  of  the  information 
reducing  aspects  of  eq.  (6),  Hfx^,  u^)  is  usually  singular  so  that  pseudo-inver¬ 
sion  must  be  used  carefully  while  avoiding  the  system  noise  usually  associated 
with  ill-conditioned  systems.  The  analogous  technique  for  space -invariant 
systems  is  to  use  Fourier  techniques  to  diagonalize  H,  thus  simplifying  the 
inversion  if  a  careful  choice  of  spectral  cutoff  is  made  to  reduce  system  noise. 
For  inversion,  singular-value  decomposition  (SVD)  techniques  [7-10lhave  been 

used  to  obtain  a  unique  pseudo-inverse  H*  (u  ,  x  )  which  has  then  used  in  the 

r  r 

restoration  operation 
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-^.(xe»uT)  =  H.+(ur,  xr)i(x0,  xr) 


(7) 


The  ^estimate  Q  (xq*  uf)  produced  by  this  procedure  minimizes  the  functional 
||  -  I  II  with  the  measured  image  I  and  known  H,  and  at  the  same  time 

A  A 

finds  the  £  with  minimum  length  ||  &  ||  .  (here  the  arguments  x  ,u  and  x 

min  °  r  r  9 

are  omitted).  This  is  accomplished  by  the  SVD  method  in  which  the  matrix 
His  expressed  as 


N 

H  =  £ 

i=l 


X2. 


u. 


T 
—  i 


(8) 


where  the  X  are  the  eigenvalues  of  H  HT(and  of  HTH),  the  u.  are  the  orthonor- 

T  “  “  “ i  _ 

mal  eigenvectors  of  H  H  ,  and  the  v.  are  the  orthonormal  eigenvectors  of  H  H. 

It  can  be  shown  that  the  pseudo  inverse  is  given  by 


where  X.,  i  =  1,  K  are  the  nonzero  eigenvalues.  The  SVD  computer  routine  of 
Golub  and  Reinsch  QO]  computes  these  in  a  numerically  stable  way;  by  a  judi¬ 
cious  choice  of  the  eigenvalue  cutoff  Xk  a  system  inversion  described  by  eq. 

(7)  is  obtained. 

This  procedure  has  been  implemented  on  the  100  X  100  point  image 
of  figure  2d  blurred  by  astigmatism.  Figure  2e  shows  the  result  of  a  polar 
coordinate  transformation  applied  to  figure  2d;  this  operation  converts  the 
degradation  to  a  space-variant  blur  in  one  direction.  Figure  2f  shows  the 
restoration  by  SVD,  in  which  the  5  eigenvalues  out  of  100  whose  magnitudes 
were  less  than  10  were  set  to  zero.  Following  restoration  of  each  line 
in  the  (x^,  Xg)  system,  an  :<  iverse  polar  coordinate  transformation  [ll  is 
used  to  produce  the  final  restoration  of  figure  2f.  Future  modifications  to 
this  method  may  use  the  SVD  for  Wiener  filtering  and  estimation  under 
various  constraints  and  noise  statistics. 

This  procedure  can  also  be  extended  for  restoration  of  the  general 
astigmatism  and  curvature  of  field  case.  A  brief  description  of  the  proce- 
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dure  involves:  (a)  a  polar  coordinati:  distortion; (b)  a  Fourier  transform  in 

the  9  direction  to  decouple  space-invariant  blur  as  a  function  of  u  as  ex- 

r 

pressed  in  eq.  (3b);(c)  a  space-variant  estimate  in  the  radial  direction  by 
p3eudo- inverse  methods  to  obtain  U (u  ,  X  )  for  each  9  spatial  frequency 
variable  X  ;(d)  an  inverse  Fourier  transform  in  the  9  variable  to  obtain  O 
(ur»ufl)‘>  and(e)  a  reverse  polar  coordinate  distortion.  Although  complicated, 
this  procedure  enables  the  general  four-dimensional  problem  of  restoration 
to  be  converted  to  a  family  of  two-dimensional  space-variant  problems  vv th 
point- spread  functions  depending  on  X. 
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4.  6  Fast-Suboptimal  Wiener  Filtering  of  Markov  Processes 
Ali  Habibi 

Wiener  filtering,  the  classical  technique  of  estimating  signals  from 
noisy  observations,  has  been  developed  primarily  for  continuous  signals. 

In  reference  [l],  Wiener  filtering  of  discrete  signals  using  various  unitary 
transformations  such  as  the  Karhunen-Loeve,  Fourier,  and  the  Hadamard 
transforms  has  been  considered  in  terms  of  computational  efficiencies.  It 
has  been  shown  that  optimal  results  can  be  obtained  with  any  unitary  trans¬ 
formation  of  the  data,  as  opposed  to  operating  on  the  data  directly,  but 
optimal  transform  domain  processing  requires  more  computations.  How¬ 
ever  if  one  is  willing  to  accept  a  small  degradation  in  performance,  the 
unitary  transformations  provide  considerable  computational  s-vvings.  The 
Wiener  filter  transformation  concept  is  extended  here  to  the  use  of  lower- 
triangular  transformations. 

Consider  an  N-dimensional  data  vector  X  =  (Xj,  ^  .  .  . ,  x^7  with  a 
covariance  matrix  £x  which  has  been  corrupted  by  a  white  noise  vector  V 
=  •  •  •  *  VN)  with  a  covariance  matrix  resulting  in  a  noisy  signal 

S  X  +  V.  Wiener  filtering  S  consists  of  premultiplying  S  by  filter  matrix 

£  -  £x+  C/1  (1) 

Using  any  linear  transformation  A  on  the  observation,  and  the  inverse  of 

the  transform  A  on  the  filtered  signal  does  not  change  the  estimated 

signal  if  the  filter  response  is  altered  accordingly,  that  is,  if  instead  of  a 

filter  matrix  G,  =  A  G  A  1  is  substituted.  If  A  is  a  unitary  operator  that 

transforms  X  to  an  uncorrelated  vector  Y  then  vector  Y  possesses  a  diagonal 

covariance  matrix.  If  A  is  unitary  it  follows  that  the  covariance  of  the 

transformed  noise  vector  W  remains  diagonal.  Thus  G  is  diagonal  and 

A  ^ 

Wiener  filtering  by  it  requires  only  N  multiplications.  However  2N" 

multiplications  are  needed  to  perform  the  transformation  of  the  input  and 
the  output  of  the  filter  by  A  and  A“*. 


Wiener  Filtering  by  Lower-Triangular  Transformation*  Given  any 
covariance  matrix  £xone  can  find  lower-triangular  matrices  L  and  L-1 
such  that  L  LT  is  diagonal  fel .  This  implies  that  Y  =  L  X  is  an  un¬ 
correlated  vector  and  transforming  by  L  and  L_1  requires  a  total  of 
2(N  /2  -  N)  multiplication  operations.  However,  using  this  transformation 

on  the  noisy  signal  S  prior  to  filtering  requires  filtering  L  S  by  G  where 

L 

-L  =  --  ^_1  (2) 

Denoting  the  covariance  matrix  of  the  signal  vector  Y  and  the  noise  vector 
-  hY  £y  and  in  the  transformed  domain,  respectively, 

-Y  =  -“X“T  (3) 

and 

-  W  =  =  av  L 1?  (4) 

Using  eqs.  (3)  and  (4)  in  eq.  (2)  gives 

-L  =  -Y(-Y  +  4  --T)_1  <5> 


Since  L  is  non-unitary,  unlike  the  Karhunen-Loeve  transformation. 

-L  18  not  diagonal.  Therefore,  the  total  number  of  opera¬ 
tions  using  the  L  transformation  requires  2[(N2/2)-N]  +  N2  operations, 
which  is  still  less  than  the  number  of  multiplications  required  using  the 
Karhunen-Loeve  transform.  However,  it  is  more  than  that  required  for 
direct  Wiener  filtering.  Table  I  summarizes  the  number  of  multiplications 
required  for  Wiener  filtering  of  a  vector  of  N  components  for  various  trans¬ 
formations.  For  Markov  processes,  the  lower-triangular  transformation 
is  used  to  design  a  suboptimal  Wiener  filter  that  reduces  the  number  of 
multiplications  significantly.  This  is  possible  since  the  L  matrix  for  a 
Markov  process  is  banded.  This  result  will  be  demonstrated  for  a  first- 

order  Markov  process;  extension  to  an  mth  order  Markov  process  is  straight 
forward. 
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When  x.,  i  1, . . . ,  N  is  a  first-  order  Markov  process,  operator  L 
is  lower-triangular,  banded  with  one  off-diagonal  band,  i.  e. , 


Transformations  by  L  and  L  1  are  accomplished  recursively  as  shown 
on  figure  1  each  using  N  multiplications.  This  also  implies  that  L  LTis  a 
banded  matrix  of  one  off-diagonal  band  which  in  turn  implies  that  G_1 
is  a  banded  matrix  of  one  off-diagonal  band. 

Filtering  L  S^by  G^is  equivalent  of  solving  the  matrix  equation 


_1 

G  Y  = 

-L- 


(LS) 


(7) 


A 

for  the  estimated  transformed  signal  Y.  An  exact  solution  to  eq.  (7)  for 
known  G^  requires  N  operations,  but  if  G^1  is  positive  definite  and  possese£ 
only  m  off-diagonal  bands  it  is  shown  that  eq.  (7)  can  be  solved  using  only 

2(m+l)N  multiplications  [3].  Matrix  G^1  can  be  approximated  by  a  positive 
definite  matrix  if  Cy  in  eq.  (5)  is  approximated  as 


— Y  S  (1_ai  (8) 

For  a  first  order  Markov  process  this  approximation  corresponds  to 
changing  only  the  first  element  of  Cy  from  unity  to  (l-a* )  which  in  turn 
corresponds  to  ignoring  the  transient  state  of  the  process.  Naturally  the 
effect  of  the  approximation  reduces  as  N  increases  or  as  a  (correlation 
of  adjacent  elements  in  the  process)  decreases.  Figure  2  shows  the  effect 
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Figure  4.  6-1.  Wiener  filtering  of  a  first  order  Markov  Process 
using  lower-triangular  transformation. 


Table  I 


Wiener  Filtering  Computation  Requirements  for 

Various  Transformations 

Transform 

Number  of 

Multiplication  Operations 

Identity 

1? - - 

Karhunen-  Loeve 

2N  +  N 

Fourier /Hadamard 

2 

N  +  2N  log  N 

Lower -Triangular 

2N2-2N 

Lower -Triangular  for  m-th 
order  Markov  "suboptimal" 

2(2m+l ) 
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of  this  approximation  on  the  theoretical  value  of  the  estimation  error  for 
various  values  of  0^  and  the  signal  to  noise  ratios  for  a  first  order  Markov 
process.  The  effect  of  this  approximation  will  be  more  pronounced  for 
higher  order  Markov  processes.  However,  the  required  number  of  com¬ 
putations  for  this  suboptimal  Wiener  filtering  is  2(2m+l)N  for  an  m**1  order 
Markov  process  which  is  even  smaller  than  the  number  of  operations  needed 
for  Kalman  filtering  of  the  same  signal. 
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4.7  Color  Image  Restoration  by  Linear  Estimation  Methods 
Clanton  E.  Mancill 

Trichromatic  color  sensing  systems,  such  as  color  film  or  a  color 
television  camera,  reduce  the  spectral  intensity  C(X)  at  an  image  point  to  a 

triple  of  numbers  (x^x^x^  through  a  set  of  three  integral  equations  of 
the  form 

X 

V  f  C(X)S  (X)dX  i  =  1, 2,  3  (1) 

\ 

The  S.(X)  functions  are  the  spectral  sensitivity  functions  of  the  three  color 
sensors.  The  colorimetric  model  of  human  color  vision  describes  the  sensa¬ 
tion  obtained  by  observing  the  color  C(X)  in  terms  of  three  numbers  (t  t  t  l 

1  2'  3* 

called  tri stimulus  values,  which  are  given  by  the  set  of  equations 
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for  i  =  1,  2,  3 


(2) 


t. 

i 


C(X)  T  (X)  dX 


The  color  matching  functions  T.(X)  describe  the  spectral  response  of  the 

human  eye  under  the  colorimetric  model.  Since  the  T  (X)  functions  can 

i 

rarely  be  expressed  as  linear  weighted  sums  of  S. (X),  the  tristimulus 

values  t.  cannot  be  determined  exactly  by  observing  the  sensor  outputs  x  . 

1  i 

Lineal1  Tristimulus  Estimation  The  purpose  of  the  work  described 
in  this  section  is  to  estimate  the  tristimulus  values  of  a  color  after  obser¬ 
ving  the  sensor  outputs.  The  approach  used  is  to  discretize  the  linear  in¬ 
tegral  eqs.  (1)  and  (2)  and  then  to  solve  the  resulting  matrix  equations  for 
an  estimate  of  the  tristimulus  values  using  several  well  known  linear  esti¬ 
mation  methods.  The  matrix  equations  corresponding  to  eqs.  (1)  and  (2)  are 

x  =  S  £  +  1  (3) 

t_  =  T  c  (4) 

where  the  3  X  n  matrices  S  and  T  are  the  discretized  sensor  and  color  match¬ 
ing  functions  (including  quadrature  formulas  for  numerical  integration  over 
n  mesh  points),  t_,  x,  and_e  are  3x1  vectors  corresponding  to  tristimulus 
values,  sensor  outputs,  and  observation  errors,  respectively,  and  c  is  the 
n  xl  vector  representing  the  color  C(X).  The  matrices  S  and  T  are  known, 
the  vectors  c_and  £  are  unknown,  the  vector  x  is  observed,  and_t  is  to  be 
estimated. 

Estimation  Methods  The  linear  estimation  methods  which  have 

been  considered  thus  far  are  enumerated  below. 

(1)  Least  Squares  Estimate:  The  rows  of  the  color  matching  matrix 

> 

T  can  be  approximated  (in  the  least  squares  sense)  by  a  weighted  sum  of  the 
rows  of  the  sensor  matrix  S.  If  the  norm  in  n- space  is  chosen  to  be  the 
Euclidean  distance  norm, 


then  the  least  squares  tristimulus  estimate  is  given  by 

A  T  T  -1 

1LS  =  IS  (SS  )  x  (6) 

This  estimate  neglects  the  observation  noise  and  is,  strictly  speaking,  a 
deterministic  solution. 

2.  Minimum  Norm  Estimate:  Another  deterministic  solution  is  pro¬ 
vided  when  the  equation  x  =  £>_c  is  solved  for  c.  under  a  linear  constraint. 
Since  the  equation  is  underdetermined,  there  are  many  solutions.  The 
solution  which  minimizes  a  norm  of  the  form 

lie  llN  =  (£TNc)2  (7) 

is  chosen,  where  N  is  a  positive  definite  matrix  chosen  so  that  ^maximizes 
some  known  a  priori  property  of  spectral  functions,  such  as  smoothness. 

For  example,  N  could  be  selected  so  that  the  solution  to  x  =  Sc  is  chosen 
which  has  the  smallest  average  squared  second  difference.  The  minimum 
norm  tri stimulus  estimate  is  given  by 

!mn=1N'1ST(SN-1ST)-Ix  (8) 

3.  Minimum  Bias  -  Minimum  Variance  Estimate:  Assume  that  the 
observation  error  _E_has  mean  and  covariance  given  by 

E(£)  =  0  E(e  eT)  =  R  (9) 

An  unbiased  linear  estimator  of_t  or  c  is  not  possible  since  the  equation 
i  =  S  c  is  underdetermined.  The  estimator  is  restricted  to  the  class  of 
minimum  bias  estimators  of  minimum  variance.  The  resulting  trisiimulus 
estimator  is  then  given  by 

4*BMV  =  I  M  ST  (S  M  sV1  X  (10) 
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rp  l 

where  the  M-norm  '■£  II  ^  =  (c  Me)1  is  used  in  the  n-dimensional  space 
and  the  Euclidean  norm  is  used  in  three  dimensional  space.  Note  that  the 
estimate  does  rot  contain  If  Mis  replaced  by  the  identity  matrix,  eq.  (10) 

becomes  eq.  (6),  the  least  squares  estimate. 

The  three  estimates  considered  thus  far  are  of  the  form 

i  IS.  -  IS+  x  (11) 

where  .S  is  a  generalized  inverse  of  S. 

4.  Wiener  Estimate:  If  the  color  vector  £  is  itself  a  random  vec¬ 
tor  whose  first  two  moments  are  known,  then  it  is  possible  to  derive  a 
Wiener  estimator  which  gives  a  solution  of  eq.  (3)  that  minimizes  the  ex¬ 
pected  squared  error 

Q  =  E[(c  -  £)T(c  -  £)  1  (12) 

The  moments  of  c^  and  the  observation  errorj  are  given  by 

E(c)  =  y.  e£(c  -  y.J(£ -yJT]  =  R.  (13) 

C  v  “  c  c 

E(e)  =0  E[e  eT]  =  R 

—  -  —  ee 

The  resulting  c,  when  multiplied  by  T,  gives  the  Wiener  tristimulus  esti¬ 
mate 


t 

— w 


T(m+R  S  (SR  S+R) 
c  —  cc- - cc  —  -ee 


-1 


x) 


(14) 


Simulation  Results  The  four  types  of  estimates  described  have  been 
tested  using  as  inputs  a  set  of  ten  test  colors.  These  represent  spectral  dis¬ 
tributions  of  ns  ural  colors  (sky,  grass,  flesh  tones,  etc,  )  seen  in  daylight. 
The  S  matrix  represents  the  layer  sensitivities  of  a  common  reversal  color 
film,  corrected  for  lens  absorption.  The  matrix  consists  of  the  color 
matching  curves  of  the  Uniform  Chromaticity  Scale  (UCS)  system.  Eighty 
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mesh  points  have  been  used  in  the  discretization  of  S(X),T(X),  and  C(X). 

The  minimum  norm  estimate  utilized  a  smoothing  matrix  N  which  mini¬ 
mizes  the  average  squared  second  difference  of  c.  The  minimum-bias  , 
minimum  variance  estimate  contains  M  =  _I  giving  an  estimate  equivalent 
to  the  least  squares  estimate.  The  Wiener  estimate  utilizes  a  first 
order  Markov  covariance  with  an  adjacent  element  correlation  of  0.  95, 

and  assumes  a  zero  observation  error.  The  mean  vector  u  has  been 

— c 

set  equal  to  the  illumination  curve  times  a  constant  reflectivity  of  0.  3. 

The  results  are  shown  in  Table  1.  The  chrominance  error  is  the  RMS 
error  distance  in  uniform  chromaticity  space  averaged  over  the  ten  test 
colors,  and  the  luminance  error  is  the  normalized  RMS  value  averaged  over 
the  ten  colors.  It  is  clear  that  specifying  more  and  more  a  priori  informa¬ 
tion  on  £  improves  the  estimate  of  £.  Further  improvement  might  be  made 
by  imposing  a  nonlinear  boundedness  constraint  on  c,  since  spectral  re¬ 
flectivities  must  lie  between  zero  and  one  at  all  wavelengths. 


Type  of  Estimate 

Chrominance  Error 
(ten  color  RMS  average) 

Luminance  Error 
(ten  color  RMS  avg. 

Least  Squares 

.  023 

.046 

Minimum  Bias -Minimum 

Variance,  M  =  I 

Minimum  Norm 

(smoothness  norm) 

.  022 

.034 

Wiener 

.  012 

.005 

Table  1. 

Linear  Tristimulus  Estimator  Luminance  and  Chrominance  Errors 
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4.8  Film  Recording  of  Color  Images  from  a  Television  Monitor 
Robert  Wallis 

The  end  result  of  computer  image  processing  is  generally  a 
photographic  image.  Unfortunately,  the  gross  distortions  of  grey  level 
and  color  information  that  are  inherent  in  film  recording  are  typically 
ignored.  The  result  is  often  an  image  which  has  been  processed  to  a 
greater  degree  by  the  film  than  the  computer.  The  goal  of  this  research 
is  the  capability  of  controlling  the  photographic  process  not  by  chemical 
means,  but  by  an  appropriate  pre-distortion  of  the  image  in  order  to 
neutralize  the  film's  distortions. 

C°lP.Ldi8Play  The  color  monitor  is  an  additive  three  color  device 
which  generates  a  spectral  distribution  given  by 

3 

=  £  P{  P.  (X) 
i  =  1  1 

where  P.(X)  i  =  1,  2,  3  are  the  spectral  distributions  of  the  red,  green,  and 
blue  phosphors  and  the  p^  i  =  1,  2,  3  are  the  weights  given  to  each  phosphor. 
If  the  nonlinearities  of  the  CRT  are  accounted  for  (gamma  correction)  the 
weights  in  the  summation  correspond  to  the  drive  signals  of  the  monitor. 
These  weights  are  often  denoted  as  R,  G,  and  B. 

Color  film  Color  film  [2l  is  difficult  to  analyze  for  two  reasons. 
First,  it  is  a  subtractive  process.  That  is,  it  generates  colors  by  blocking 
certain  wavelengths  from  white  light,  whereas  an  additive  system  super¬ 
imposes  various  spectra.  The  dyes  used  as  subtractive  primaries  are 
cyan,  magenta,  and  yellow.  These  can  be  considered  notch  filters  which 
reject  red,  green  and  blue  respectively.  The  use  of  red,  green,  and  blue 
primaries  for  a  subtractive  system  would  be  unfeasible.  This  is  because 
red,  green,  and  blue  are  generated  by  passing  narrow  spectral  portions  of 
white  light.  Thus,  if  any  two  such  dyes  were  to  be  mixed,  the  resultant  dye 
would  block  nearly  all  wavelengths,  since  their  passbands  would  not  overlap. 

A  second  difficulty  is  the  nonlinear  response  of  the  photographic 
emulsion  to  light.  The  relationship  between  exposure  and  the  resulting 
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optical  density  is  known  as  the  Hurter-Driffeld  curve  fe]  .  If  the  image 
intensity  at  each  point  of  the  image  to  be  photographed  is  predistorted  by 
the  inverse  of  this  curve,  the  film  can  be  forced  to  yield  a  linear  response. 

Tr.istimulus  values  The  basis  of  colorimetry  is  the  three  color 
theory  of  human  vision,  which  states  that  human  perception  of  color  is  in¬ 
trinsically  a  three  dimensional  phenomenon.  Specifically,  for  each  spec¬ 
tral  distribution  S(X)  color  perception  can  be  characterized  by  three  tristi¬ 
mulus  values  given  by 

t.  =  J*T.(X)S(X)dX  i  =  1,2,3 

where  the  T  (»  are  known  as  color  matching  functions.  The  t  are  scaled 

i 

such  that  t.  =  1  for  a  reference  white.  In  order  for  two  colors  to  match, 
it  is  sufficient  that  their  tristimulus  values  match,  but  it  is  unnecessary 
that  their  spectral  distributions  match.  This  of  course  presupposes  that 
the  colors  are  compared  under  similar  circumstances.  Tristimulus  values 
can  be  transformed  into  a  number  of  different  coordinate  systems  in  order 
to  facilitate  interpretation.  One  system  which  has  been  found  useful  is  the 
cube  root  coordinate  system  [3],  which  is  described  in  figure  1.  The  cube 
root  color  space  has  the  advantage  that  equal  distances  in  the  spuce  very 
nearly  correspond  to  equal  distances  perceptually.  Thus,  it  can  be  utilized 
to  judge  colorimetric  errors  between  pairs  of  colors. 

Computer  simulation  of  color  film  Using  the  model  of  figure  2,  a 
simulation  of  color  transparency  film  has  been  performed.  The  hypothetical 
inputs  are  the  three  individual  color  monitor  phosphors  exposed  over  a  wide 
range  of  "f- stops".  (Each  f-stop  corresponds  to  a  doubling  of  the  exposure.  ) 
The  results  are  given  in  figure  3.  The  loops  are  the  paths  in  color  space  that 
retult  as  the  exposure  is  incremented  in  one-half  f-stop  intervals.  In  an 
attempt  to  verify  the  accuracy  of  the  model,  red,  green  and  blue  fields  of 
a  Conrac  color  monitor  have  been  photographed  using  Ektachrome-X  film 
over  a  similar  range  of  exposures.  The  resulting  colors  of  the  transparencies 
have  been  measured  using  a  colorimeter.  The  agreement  between  theory  and 
experiment  is  good,  as  can  be  judged  from  figure  3.  It  is  expected  that  more 
careful  modelling  will  yield  even  better  correlation.  Evidently,  the  exposure 
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L  =25.29  [(I00Y)I/3-  18.33] 
a  =106.0  [( 102 X)l/3  —  ( I00Y  )l/3] 

b  =  42.34  [OOOY)l/3-(84.7Z)'/3] 

where  X,Y,  Z  are  the  C.I. E.  tristimulus  values 


Figure  4.8-1.  Cube  root  color  coordinate  system. 
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is  a  parameter  of  great  concern,  in  that  large  shifts  in  hue  and  saturation 
result  from  even  a  one-half  f-stop  change  in  exposure.  For  instance,  the 
red  shifts  towards  orange  an/  the  blue  shifts  towards  cyan  as  exposure  is 
increased.  Another  general  result  is  the  loss  of  saturation  as  compared  with 
the  original  color  photographed.  Note  the  greater  saturation  of  the  original 
red,  green  and  blue  primaries  of  figure  3.  A  simple  algorithm  providing 
a  first  order  correction  to  this  loss  of  saturation  has  been  developed.  It 
essentially  subtracts  white  (or  neutral  density)  from  each  color  to  be  photo¬ 
graphed,  but  maintains  the  original  luminance.  Thus,  the  loss  of  saturation 
suffered  by  film  recording  is  somewhat  offset  by  the  pre- saturation.  The 
algorithm  is 


Y  =  .  241  R  +  .682  G  +  .  077  B 
a  =  k  Y 

where  Y  is  the  luminance  of  the  color  and  k  is  a  constant  which  determines 
the  degree  of  presaturation  desired.  The  R,  G  and  B  are  then  transformed 
as  follows 


R'  ■  <t^t>  (R  •  •* 

G'  ■  ^ri(G  -  <* 

B’  =  <t (B  -  °> 

Note  that  the  luminance  remains  invariant,  but  that  the  differences  between 
R,  G  and  B  are  exaggerated.  The  algorithm  has  been  tested  on  an  actual 
color  image,  and  has  been  found  to  result  in  a  marked  improvement.  An 
algorithm  providing  a  colorimetrically  exact  correction  is  possible,  assuming 
the  desired  color  is  within  the  film's  gamut  of  color  reproducibility.  That  is, 
given  a  desired  output  color,  the  required  input  color  can  be  found  exactly 
by  iterative  techniques.  However,  for  a  large  image,  the  required  compu¬ 
tation  time  would  be  prohibitive.  A  compromise  algorithm  employing  look 
up  tables  is  being  developed  which  will  combine  speed  and  accuracy. 
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4.  9  Psycho- Visual  Pseudo- Color  Mapping  of  Imagery 
Warner  Frei  and  Mark  Stein 

Pseudo- color  techniques  generally  map  two-dimensional  functions 
of  space  variables  into  a  pre-selected  gamut  of  colors,  thus  creating  arti¬ 
ficially  or  pseudo- colored  images  for  easier  and  more  efficient  information 
extraction  by  human  observers.  These  techniques  can  be  broadly  classified 
into  two  categories:  (a)  mappings  of  a  scalar  spatial  function  such  as  low- 
contrast  achromatic  images  into  some  fixed  arbitrary  color  scale  so  that 
minute  differences  of  grey  in  the  original  image  may  then  appear  as  easily 
recognizable  color  differences  &l;(b)  mapping  of  a  vector-valued  spatial 
function  into  color  space.  Since  the  human  visual  system  is  able  to  recog¬ 
nize  three  attributes  of  colors  mon  or  less  independently  of  each  other 
(attributes  such  as  brightness,  hue  a*?d  saturation)  one  can  theoretically 
map  three  components  of  a  vector-valued  spatial  function  into  three  suitably 
chosen  attributes  and  instruct  the  observer  that  for  example,  brightness 
corresponds  to  function  A,  h  ie  to  function  B  and  saturation  to  function  C. 

It  is  felt  that  many  pseudo-color  techniques  often  fail  to  re«ch 
their  potential  because  they  are  designed  with  little  understanding  of  the 
visual  system.  Arbitrary  color  assignments,  indeed,  tend  to  introduce 
artificial  contours  and  destroy  other  valuable  features  creating  more  confu¬ 
sion  than  clarification.  Careless  pseudo-coloring  may  thus  result  in  pictures 
with  a  mere  artistic  interest.  In  this  report,  these  points  will  be  illustrated 
by  an  application  of  pseudo- color  to  nuclear  medicine  imagery. 

Comparison  of  Inhalation  and  Perfusion  Lung  Scintigrams  Lung 
scintigrams  or  scans  are  images  of  the  lung  obtained  by  imaging  the  output 
of  a  Y-  ray  detector  which  is  placed  over  a  patient  who  has  received  an  in- 
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ternally  administered  radio  nuclide.  In  lung  scanning,  two  images  are  obtained 
which  reflect  the  lung  blood  pervasion  and  ventilation  respectively.  An  area 
that  is  normally  ventilated  but  not  perfused  most  likely  has  blood  clot  in  the 
vessels  (pulmonary  embolus).  Tnirt  distinction  is  extremely  important  in 
clinical  therapeutic  decision-making.  Since  the  V  -  ray  camera  is  an  imaging 
system  with  a  low  spatial  resolution  and  the  images  are  extremely  noisy,  the 
correspondence  determination  between  the  two  images  is  very  difficult  for 
human  observers,  especially  when  multiple  defects  are  present  on  both  types 
of  scintigrams.  Identification  of  this  correspondence  requires  the  ability  to 
determine:  (a)  degree  (severity)  of  defects;  (b)  polarity  of  defects  (which 
image  has  the  defect;(c)  anatomic  limits  (lung  outlines).  Viewing  the  per¬ 
fusion  and  inhalation  scans  adjacent  to,  or  superimposed  upon  one  another 
results  in  an  inability  to  extract  these  three  information  elements  accurately. 

The  perfusion  and  inhalation  scans,  denoted  as  two  functions  f  (x.,  x2) 
and  fj(Xj,  x^)  of  the  same  space  variables  x^,  x^,  contain  combined  anatomic 
(structural)  and  defect  information.  In  order  to  combine  the  information  of 
interest  into  one  color  image,  a  field 


Y(x1»x2)  =  kj  max  ^fp(xi'x2)'fi(xi»x2)  k2 


(1) 


which  preserves  structural  information,  and  a  field 


(x1»x2)  =  M p(xl,x2)  -  V’VXjjH 


(2) 


which  measures  differences  in  count  rates  are  formed  (k  ,  k  ,  k  are  arbi- 

X  C*  J 

trary  scaling  factors).  Next,Y(Xj,  x^)  is  mapped  into  luminance  and  6(Xj,x2) 
into  an  axis  of  complementary  colors.  In  the  general  case,  such  an  axis 
is  defined  in  the  CIE  chromaticity  diagram  by  a  straight  line 


y  =  ax  +  b 


(3) 


y  =  ax  +  b 
w  w 


Choosing  a  zero  difference  (normal  area)  to  be  represented  by  white. 
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(*w  -  Yw  =  1/3),  the  chromaticity  of  the  displayed  pixel  is 


x(6)  = 

1  /  3  +  6  (1  +  a2)  '2 

(4) 

y(6)  = 

-2  -  - 

1  /  3  +  6  (1  +  a  )  2 

(5) 

Since  the  chromaticity  coordinates  x,  y  are  given  by 

=  X  _  Y 

X  X+Y+Z  y~X+Y+Z 

where  X,  Y,  Z  are  tristimuli  referred  to  the  CIE  non-physical  primaries 
(X),  (Y),  (Z),  one  obtains 


Y 

X 

Z 


kj  max  tfp(x1#x2),  fjfxj,  x2)]+  k2 

Y  *ii 
y(6) 

1  -  x(6)  -  y(6) 

y(6) 


(6) 

(7) 


where  x(£),  y(6)  are  given  by  eqs.  (4)  and  (5).  Finally  the  signals  to  drive  a 
color  TV  monitor  are  obtained  by  a  change  of  coordinates  in  tri stimulus 
space  as  given  by 


’r 

= 

'ail 

ai2 

a13* 

G 

= 

a21 

a22 

a23 

B 

»  « 

= 

a31 

• 

a32 

a33 

p  • 

X 

Y 

Z 


(8) 


A  much  more  computationally  economic  solution  is  obtained  by  restricting  the 
choice  of  complementary  colors  to  an  axis  passing  through  one  of  the  receiver 
primaries.  Considering  further  the  energy  limitations  of  the  display  pri¬ 
maries,  it  is  found  that  the  most  efficieit  axis  is  the  green-magenta  axis, 
because  the  luminous  efficiencies  of  the  NTSC  primaries  give  y  =0.  585 
and  yR  +  yB  =  yM  =  °*  413.  The  pseudo-color  mapping,  such  that  6  =  0 
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maps  into  white  becomes 


and 


G 

G  +  M 


6+  0.  5 


Y  =  0.  585G  +  0.413  M 


(9) 

(10) 


whure  G(xJ,x2)  and  Mfx^x^  =  Rfx^x^  =B(x1,x2)  are  the  green,  red  and  blue 
signals  for  a  TV  display.  Letting 


K  =  M—J  for  -0.  5  <  6  ^  0.  5 

0.  5  +  o 


gives 

G(x  ,  x2)  = 

YfXj,  x2) 

0.  585  +  0.413 

(ID 

and 

M(xlf  x2)  = 

kG  -  R(Xj,  x2)  =  B(Xj,x2) 

(12) 

With  the  mapping  described,  pseudo-color  images  are  obtained  where  inten¬ 
sity  naturally  carries  outline  information,  hue  indicates  the  polarity  of  the 
defects,  and  saturation  shows  severity,  while  normal  areas  appear  white. 

The  constant  k2  insures  that  the  intensity  does  not  drop  below  a  level  that 
would  make  color  discrimination  difficult. 

Initial  studies  have  shown  that  the  ability  of  the  physician  to  deter¬ 
mine  accurately  correspondences  in  perfusion  and  ventilation  scans  can  be 
greatly  enhanced  by  this  single  pseudo-colored  images  as  compared  to  two 
separate  black-and-white  images.  In  addition  this  rational  use  of  color 
may  increase  the  detectability  of  defects.  Both  results  have  been  recognized 
in  the  initial  work,  but  more  experience  and  rigorous  psycho -physical  testing 
will  be  needed  to  develop  confidence  of  the  medical  community  in  this  method. 

References 
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5.  Image  Data  extraction  Projects 

Image  data  extraction  is  a  name  given  to  a  collection  of  projects  that 
are  concerned  with  the  detection  of  features  within  an  image  and  methods 
of  measuring  these  features. 

A  hybrid  optical-digital  feature  extraction  system  for  images  is 
described  in  the  first  report.  The  system  includes  a  two  dimensional 
Fourier  transformation  stage  implemented  by  optics  followed  by  a  diffraction 
pattern  electronic  sensor  which  supplies  digital  signals  to  a  computer 
processor.  The  hybrid  system  can  be  programmed  for  edge  detection, 
texture  detection,  or  to  detect  a  variety  of  other  image  features. 

In  the  next  report  a  description  is  given  of  an  interactive  edge  trac¬ 
ing  device  for  photographic  transparencies.  The  device  consists  of  an 
electro-acoustic  sensor  that  provides  coordinate  points  on  an  edge  contour 
which  are  then  formatted  for  remote  computer  entry. 

The  following  two  reports  are  concerned  with  techniques  for  obtaining 
two  and  three  dimensional  perspective  views  of  an  objer'  from  one  or  two 
dimensional  density  profiles.  Methods  based  upon  the  Fourier  transform, 
convolution,  algebraic  reconstruction,  and  Kalman  filtering  are  explored. 

In  the  last  report  a  method  of  image  boundary  estimation  is  discussed. 
Experimental  results  are  given  for  several  test  images. 

5.  1  An  Optical-Digital  System  for  Feature  Extraction  From  Images 

Richard  P.  Kruger  and  Ernest  L.  Hall 

A  hybrid  approach  to  image  feature  extraction  is  described  which 
combines  the  parallel  processing  capability  of  a  coherent  optical  system 
with  the  logical  processing  ability  of  a  computer  controlled  scanner.  A 
novel  feature  of  the  system  is  provided  through  the  use  of  an  automated 
film  transport  which  permits  parallel  combinations  of  optical  and  digital 
processing  on  local  or  global  areas  of  an  image  rather  than  a  cascade  of 
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sequential  optical-digital  or  digital-optical  operations.  This  system 
configuration  lias  been  motivated  by  the  desire  to  combine  both  digital 
and  optical  processing  capabilities  in  a  single  system,  thus  permitting 
a  convenient  switch  from  one  system  to  the  other  when  desired  for  a 
particular  task.  For  example,  the  digital  system  may  be  used  to  locate 
objects  and  the  optical  system  to  measure  properties  of  the  objects. 

A  block  diagram  of  the  hybrid  system  is  shown  in  Figure  1.  At 
present  a  film  transparency  is  loaded  into  a  fixed  frame  which  may  be 
manually  registered  in  both  the  optical  and  digital  systems.  A  future 
modification  will  permit  access  to  the  film  via  the  film  transport  shown 
in  Figure  2a.  The  film  scanner  is  an  image  dissector  type  as  shown  in 
Figure  2b.  The  principal  components  of  the  optical  system  are  the 
sampling  sensors  shown  in  Figure  3a  and  the  electronic  amplifier  unit 
shown  in  Figure  3b.  The  controlling  computer  is  a  PDP  11/10  which 
is  interfaced  via  an  HP  2100A  to  an  IBM  360/44. 

For  a  practical  application  of  the  hybrid  system,  the  transformations 
and  operations  should  use  the  particular  device  to  advantage.  For  example, 
Fourier  transformations  are  easily  computed  with  the  coherent  optical 
system.  Nonlinear  operations  such  as  threshold  computations  are  easily 
performed  with  the  computer  controlled  scanner.  Furthermore,  for  a 
feature  extraction  system,  a  large  computer  storage  is  not  required  since 
the  film  is  an  excellent  read  only  memory  and  the  number  of  required 
measurements  is  usually  small. 

In  the  hybrid  system  a  diffraction  pattern  sampler  is  used  for  the 
optical  measurements.  The  Fourier  transform  of  f(x,  y)  is  defined  as 

00 

F(u,  v)  =  f(x,  y)exp{-(uv+vy)}dxdy 
and  the  k**1  measurement,  m  ,  is  given  by 

K 

rru  =  JJ  I  F(u,  v)  |  dudk 
Rk(u,v) 
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Figure  5.  1-1.  Block  diagram  of  hybrid  optical /digital  image  processing  system. 


(b)  image  digitizer. 


(Figure  5.  1-3.  Optical  transport  and  digitizer. 
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where  the  regions  Rk(u,v)  are  annular  ring  or  wedge  shaped  areas  in  the 
spatial  frequency  domain  [l]  .  These  measurements  have  been  found 
useful  in  many  applications.  By  the  use  of  the  film  transport,  any  desired 
region  of  the  image  may  be  sampled  in  the  above  manner. 

The  digital  computation  is  completely  programmable  and  limited 
only  by  the  resolution  of  the  scanning  device  and  speed  of  the  computer. 
Typical  operations  might  include  a  low  resolution  scan  to  produce  a  computer 
image,  f(i,  j),  computation  of  the  histogram,  h(f)  and  spatial  profiles 

*(*>  = 

j 

«P(y)  =  £f(i,j) 
i 

which  are  used  to  locate  a  center  point  (x,  y)  for  the  optical  measurements. 

Another  set  of  digital  measurement  which  are  presently  used  are 

spatial  moments  [2]  m  (f)  defined  bv 

pq  7 

00 

mpq(f)  =  xPx<lf(x*  y)dxdv 

The  moment  computation  may  be  normalized  to  obtain  translation  or  size 
invariance.  The  moments  may  also  be  extended  to  reflect  edge  structure 
of  texture  patterns.  F or  example,  the  picture  function  at  any  point  (?,  T]) 
may  be  expressed  as  a  Taylor's  series  about  the  point  (x,y).  The  first 
three  terms  of  this  expansion  are 

f(x+?,  y+n)  =  f(x,y)  +  g3i(xi  y)  +  n 

9x  dy 

Therefore,  the  picture  function  at  the  point  (?,  71)  may  be  expressed  in  terms 
of  the  moments 

”pqfc+?.y+H)l  =  mpq[f(x.yn+Smpqryx,yn+Hm  [£  (x,y)] 

Tne  operations  which  combine  the  optical  and  digital  measurements 

appear  to  be  highly  application  dependent  and  will  be  studied  for  particular 
applications. 
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A  current  application  of  this  hybrid  system  involves  a  hybrid 
system  for  automatic  screening  of  chest  radiographs  to  detect  the 
presence  or  absence  of  coal  workers  black  lung  disease  (pneumoconiosis). 

The  proposed  system  incorporates  coarse  digitization  of  the  entire  chest 
x-ray  film  to  either  automatically  reject  a  film  with  improper  exposure 
or  detect  and  measure  gross  anatomical  features  in  the  film.  In  addition 
aperture  centers  in  all  six  lung  zones  are  computed  for  transmission  to 
the  computer  controlled  film  transport. 

Figure  4  depicts  the  output  of  this  first  digital  process.  Shown  are 
least  squares  oolynomial  estimates  of  the  lung  outlines,  rectangular  estimates 
of  the  six  lung  zones,  and  simulated  aperture  centers  for  each  zone.  The 
second  stage  of  processing  will  involve  the  movement  of  the  film  transport 
to  each  of  the  aperture  centers,  measurement  of  their  respective  spectral 
content  and  conservative  initial  diagnostic  classification.  If  the  film  is 
classified  normal  in  all  zones,  it  will  be  removed  from  the  system.  If 
this  conservative  initial  classification  indicates  a  possible  abnormal  situation, 
the  random  access  capability  of  the  digital  scanner  will  be  used  to  raster 
scan  the  suspicous  rectangular  zones  at  high  resolution.  This  third  sta3e 
of  processing  will  involve  the  computation  of  textural  measures  obtained 
within  the  boundaries  of  the  interpolated  least  square  estimates  of  lung 
zones  in  question.  These  textural  measures  will  be  used  to  obtain  a 
second  diagnostic  classification.  This  final  decision  will  either  place  the 
film  into  the  normal  group  or  place  it  in  one  or  se  /eral  disease  classes. 

This  flexible  hybrid  approach  allows  for  processing  economy  at  each 
stage  of  analysis.  In  general  use  of  the  film  as  a  read  only  memory,  different 
optical  detectors,  random  access  digital  scanning  and  a  computer  controlled 

film  transport  will  provide  a  flexible  laboratory  tool  for  image  analysis  and 
data  extraction. 
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5.2  A  Terminal/Timeshare  Based  Device  for  Interactive  Boundary 
Tracing  and  Analysis 
Richard  P.  Kruger 

The  ability  to  input  graphic  information  such  as  image  boundary 
information  and  strip  chart  data  in  digital  form  into  a  computer  for  analysis 
can  be  approached  in  mar.y  ways.  The  specific  facility  to  be  described  here  is 
particularly  applicable  in  circumstances  whi  h  preclude  the  use  of  a  local 
minicomputer  system,  and  where  it  is  desired  to  share  a  common  program 
library  with  several  remote  users. 

A  system  with  these  general  attributes  has  been  installed  at  the  USC- 
LA  County  Medical  Center  and  connected  to  the  USC  main  campus  computer 
via  standard  dial  up  telephone  lines.  The  library  consists  of  Fortran  language 
programs  accessed  under  IBM  time  share  option  for  analysis  of  various 
aspects  of  left  ventricular  heart  function.  The  system  configuration  for 
data  collection  is  shown  in  Figure  1.  The  digital  input  device  consists  of 
a  graf  pen  digitizing  tablet. 

The  Graf-Pen  graphic  tablet  digitizer  utilizes  acoustic  ranging  to 
digitally  measure  the  position  of  a  stylus  on  the  surface  of  the  tablet.  The 
measurement  is  performed  by  measuring  the  time  required  for  sound  to 
travel  from  a  spark  sound  source  on  the  tip  of  the  stylus  to  a  pair  of  linear 
microphones  located  along  two  axes  of  the  tablet.  This  time  measurement 
is  translated  into  a  number  between  zero  and  1999  for  each  axis,  effectively 
dividing  the  tablet  into  a  matrix  of  2000  x  2000  points. 

The  four  digit  binary-coded  decimal  numbers  generated  by  the  Graf- 
Pen  for  each  data  point  are  converted  to  a  sequence  of  characters  suitable 
for  transmission  to  a  computer  terminal  or  teletype  unit  by  a  data  converter. 

1  his  data  converter,  in  addition  to  performing  the  parallel  to  serial  data 
conversion,  inserts  control  characters  into  the  data  string  to  make  it 
compatable  with  a  computer. 
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Figure 


SERIAL  DATA  STREAM 
WITH  CONTROL  CHARACTERS 


5.2-1.  Data  acquisition  and  processing  system. 
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In  the  present  installation,  data  generated  by  the  Graf  Pen  is  passed 
from  the  data  converter  to  a  Novar  computer  terminal  which  is  equipped  with 
a  cartridge  data  tape  recorder.  This  tape  recorder  is  used  as  a  data  buffer 
to  hold  pen  data  until  it  is  required  by  a  computer  program  and  to  provide 
for  maximally  efficient  transmission  of  data  to  the  computer  during  program 
execution.  The  sporadic  manner  in  which  time- shared  computers  schedule 
program  execution  makes  such  a  buffer  a  necessity.  During  the  execution 
the  operation  of  the  tape  recorder  is  controlled  automatically  by  the  computer 
making  manual  intervention  unnecessary. 

Data  transfer  rates  for  the  graphic  data  terminal  are  determined  by 
the  rate  limitations  of  the  computer  terminal  and  the  computer's  input  equip¬ 
ment.  Since  the  existing  equipment  is  limited  to  no  more  than  30  characters 
per  second,  the  data  acquisition  and  transfer  rate  is  limited  to  approximately 
two  X-Y  coordinate  pairs  per  second.  During  digitization  operations  the 
digitizer  can  be  set  to  acquire  data  only  when  the  operator  indicates  a 
desired  point  by  pressing  on  the  data  stylus. 

In  operation  a  stop  motion  movie  projector  image  of  the  dye  filled 
left  ventricle  is  projected  onto  the  graf  pen  surface  where  the  operator 
digitizes  30  to  50  perimeter  points  obtained  from  2  to  60  frames  of 
cine-film  per  patient.  This  boundary  information  as  well  as  digitized 
strip  chart  recording  of  a  simultaneously  recorded  pressure  versus  time 
curve  is  stored  on  cassette  tape  for  transmission. 

At  this  point  the  IBM  time  share  option  is  accessed  with  resident 
programs  to  accept,  analyze  the  data  and  return  a  report  to  the  operator. 

The  programs  initially  interact  with  the  operator  in  a  question  and  answer 
mode  followed  by  computation  and  return  of  20  quantitative  indices  of 
ventricular  function.  In  addition,  a  high  resolution  (30  characters  per  inch) 
plot  routine  returns  a  plot  of  the  end  diastole  (full  expansion  and  end  systole 
(full  contraction)  boundary  trace  (see  Figure  2)  as  well  as  volume  vs  time, 
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pressure  vs  volume,  and  force  vs  length  plots  where  applicable.  This 
information  then  becomes  part  of  the  permanent  patient  record. 

This  interactive  terminal  based  device  has  significantly  reduced  the 
time  to  complete  such  an  analysis  to  the  point  of  practical  cost-effectiveness 
The  modest  cost  of  such  a  system  will  allow  many  remote  users  to  have  ac¬ 
cess  to  a  large  general  purpose  computer  where  they  can  use  as  well  as 
augment  existing  ventricular  analysis  programs  for  use  by  the  user 
community.  This  concept  has  the  potential  for  encouraging  resource 
sharing  within  its  limited  scope. 

5.  3  Reconstruction  of  Three-Dimensional  Object  Arrays  From  Their 

Projections 

Alexander  A.  Sawchuk  and  Ernest  L.  Hall 

A  problem  of  some  interest  in  image  processing  concerns  the  re¬ 
construction  of  arrays  of  data  from  projection  measurements  having  a  lower 
spatial  dimensionality.  Typically,  the  objective  is  to  reconstruct  three- 
dimensional  object  arrays  from  two-dimensional  projections,  or  to 
reconstruct  two-dimensional  object  arrays  from  one-dimensional  projec¬ 
tions.  The  technique  has  applications  in  a  wide  variety  of  areas,  not  all 
of  which  are  traditionally  connected  with  image  processing.  Some  of  these 
application  areas  include:  radiography,  or  X-ray  imaging  [l,2];  nuclear 
scintillation  imaging  [3];  electron  microscopy  [4];  electromagnetic  ana 
pressure  wave  field  mapping  such  as  in  radar,  radioastronomy  and  sonar  [5^] 
measurement  of  volume  refractive  index  and  temperature  profiles  [6];  and 
multidimensional  picture  representation  and  coding  [7]. 

A  large  number  of  different  techniques  have  been  proposed  for 
accomplishing  image  reconstruction.  In  its  most  general  form,  the  M- 
dimensional  projection  vector  £  is  obtained  by  a  linear  combination  of  the 
N -dimensional  orginal  object  vector  f_  as  expressed  by 

£  =  Hi  (1) 
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where  H  is  an  M  XN  matrix.  All  reconstruction  algorithms  attempt  to 
produce  an  estin.ate  £ _of  the  true  There  are  many  obvious  similarities 

to  the  image  restoration  problem,  and  the  special  nature  of  reconstruction 
has  led  to  several  general  approaches:  Fourier  transform  [l,4]  ,  convolution 
[5,  8],  and  algebraic  [2].  No  one  approach  has  been  found  superior  to  the 
others,  and  extensive  analytical  comparison  is  difficult.  Other  problems, 
such  as  measurement  noise  in  taking  the  data  according  to  eq.  (1),  and 
imperfect  knowledge  of  H  also  led  to  different  tradeoffs  between  techniques. 
Noise  is  utlimately  the  most  significant  problem  in  image  processing,  and 
very  few  reconstruction  methods  really  consider  it  adequately. 

The  Fourier  transform  method  of  reconstructions  is  one  of  the  simplest 
techniques  although  it  requires  that  projection  data  be  taken  in  a  full  180° 
around  the  object  to  be  imaged.  The  method  relies  on  the  fact  that  the  Fourier 
transform  of  a  two-dimensional  projection  of  a  three-dimensional  trans¬ 
form  of  the  object.  If  ffrj,  x2>  x^  represents  the  object,  the  three-dimensional 
Fourier  transform  is  given  by 


— (fl'i2,f3)  =  JB*  ^(xl'x2'X3)expt_2Vflxl  +  f2X2  +  f3*3)ld*1d*2dx3  (2> 

and  the  central  section  of  the  transform  is 
00 

— (fj.'  *2’  0)  =  JJJl(xi*x2*x3)exp[-2TTj(f1x1  +  f2x2)]dx1dx2dx3  3) 

After  performing  the  integration  with  respect  to  x3  in  eq.  (3)  and  identifying 

00 

—3  (xl '  X2}  =  (»-(xrX2'X3)dx3  (4) 


as  the  projection  on  the  Xj-x2  axes  and  Fff^.O)  as  its  transform,  the  full 
three-dimensional  transform  can  be  built  up  plane  by  plane  using  the  trans¬ 
forms  of  different  projection  views  oi  the  object.  It  can  be  shown  that  a 
coordinate  system  (Sj,  s2)  rotated,  as  t>hown  in  figure  la,  by  6  with  respect 
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to  (x1#x2)  moves  in  exactly  the  same  manner  as  the  rotated  Fourier 
transform  space  Thus,  projections  are  taken  at  different  6 

orientations,  transformed,  and  inserted  into  the  two-dimensional  transform 
plane  shown  in  figure  lb.  The  entire  transform  of  the  object  can  be  built 
up  plane -by-plane  by  this  method.  Since  only  a  finite  number  of  projections 
are  available,  some  interpolation  in  the  transform  space  is  required. 
Several  methods  including  pulse  approximation,  linear  interpolation, 

Fourier  series  interpolation,  and  sine  function  expansion  for  band  and 
space-limited  objects  are  being  considered,  and  tradeoffs  and  relative 
accuracies  are  under  study.  Related  to  this  are  problems  of  appropriate 
sampling  rate  and  projection  density  for  adequate  reconstruction  as  a 
function  of  object  parameters.  Although  the  Fourier  transform  method 
needs  180  of  data,  it  is  relatively  easy  to  implement  and  gives  some  of 
the  best  quality  results  with  relatively  low  required  computation  when  the 
projection  data  is  not  too  noisy. 

Tne  second  major  reconstruction  technique  is  known  as  the  convolu¬ 
tion  method  [5,  8J  and  operates  entirely  in  the  spatial  domain  without  the 
use  of  any  transforms.  Denoting  the  desired  object  to  be  reconstructed  in 

polar  coordinates  by  f(r,  9),  then  a  forward  Fourier  transform  relationship 
to  transform  F(R,  $)  is 

F(R,<)  =  J02Yf(r-e)e'j2nRCO8,e',)^9  (5) 

with  reverse  transform 

Hr,  9)  =  A"  F(R,»)e+j2mRco8(9-*Wd«.  ,6) 

Equation  (6)  can  also  be  expressed  as 

A 

f(r.9)  =  p’2  |R|F(R,})e+j2T,rRcos'e-«dRd»  <7| 

"T 
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where  the  finite  limits  on  R  represent  the  bandlimited  nature  of 

F(R,  $).  Equation  (7)  is  the  basis  for  the  terhnique. 

The  method  begins  by  first  measuring  projections  g(£,  $)  as  a  function 

of  the  variable  l  at  angle  $as  shown  in  figure  2a.  This  can  be  expressed  as 
an  inverse  transform 

A 

g(M)  =  J  F(R,$)e+j2TTR*  dR  (8) 

A 
“  2 


in  a  frequency  space  identical  to  that  shown  in  figure  lb.  The  second  step 

is  to  linearly  filter  g(Z,  I)  with  a  space -invariant  impulse  response  h$(i) 
given  by 


h,«> 


/  i.i 

A 
"  2 


j2lTR£ 
eJ  dR 


(9) 


to  produce  g'(£,  $)  according  to  the  convolution  relationship 


g'(i,$)  =  g(*,§)*hf(*)  (10) 

This  operation  is  known  as  "rho-filtering"  due  to  the  radial  spatial  frequency 
amplification  of  eq.  (9)  shown  in  figure  2b.  By  the  convolution  theory  of 
transform  theory,  eq.  (10)  is  the  same  as 


A 

=  iA  |R|  F(R,$)ej2TTR<!dR  . 

"T 

The  final  step  is  to  use  g'(£,  i)  and  compute 

TT 

Ur,  9)  =  J  g'(r  cos(9-$),  $)d$ 

0 


(11) 


(12) 


-111- 


mjm 


to  give  the  object  density  f(r,  9).  Substituting  eq.  (11)  into  eq.  (12)  gives 
the  identity  of  eq.  (7)  and  proves  the  method.  In  a  discrete  sampled -data 
format  the  operations  of  eq.  (10)  and  eq.  (12)  becomes 

00 

g’(na,  9/  =  a  ^  g(ma,  9)h((m-n)a)  (13) 

m=-“ 


and 


N 

f(r,§)  =  f(jr0,k$0)  =  ^  g'[jrQ  cos(k§0-t60),teo]  (14) 

respectively.  Other  alternative  techniques  may  involve  reversing  the  order 
of  the  operations  of  eqs.  (10)  and  (12),  or  combining  them  with  part  of  the 
data  taking  operation  of  eq.  (8).  As  with  the  transform  method,  the  discrete 
nature  of  the  computer  implementation  requires  interpolation,  sampling,  and 
numerical  integration;  work  is  underway  to  study  the  best  method  of  performing 
these  tasks.  Some  study  of  the  effects  of  noise  is  also  being  considered  for 
this  system,  and  improvements  to  the  inverse  filtering  operation  in  the  rho 
filter  are  possible. 

The  third  reconstruction  technique  may  be  called  the  algebraic  approach. 
The  problem  is  illustrated  in  figure  3  and  may  be  formulated  in  the  following 
manner.  Define  a  linear  set  of  projection  equations 

Hi  =  £ 

where 

T 

i  =  fj^)  N  =  n  X  n 

T 

£  =  (gr  g2,  •  •  • ,  gN,  8N+2»  •  •  •  .  g2N»  •  •  •  »  82N+2’  g2N’  *  *  ‘  gK) 

with  K  =  M  -  (m  -1)  linearly  Independent  projection  equations,  and  H  is  a  K  x  N 
matrix  depending  only  upon  the  geometry.  These  equations  are  exactly  similar 
in  form  to  the  restoration  equation  and  only  differ  by  construction.  Some  of 
the  possible  solution  constraints  include: 
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RECTANGULAR  OR  FAN  BEAM 
GEOMETRY  ma  ANGLES 


Ray  j 
at  angle  6 


mm 


a)  Minimum  Energy 


Ji<.fJ  =  ltl 

b)  Smoothness 

J  (f)  =  fTCf 

Z  "  “*  — 

c)  Minimum  mean  square  error.  Given  a  measurement  with  noise 


1  =  &  +  i 

J30  =  111- Hill 


Solutions  may  easily  be  developed  for  these  minimization  problems. 
F or  example,  one  may  minimize 


J(f)  =  fTf  +aiTCf 


subject  to  the  constraint 


Hf  =  £ 


which  gives 


L  =  -2[H(I  +  O.C)’1HT]"1  £ 


Similarly,  one  may  minimize 


J(f_)  = 


subject  to  the  noise  measurement  constraint  which  gives 

1=  [l  +  hth  +  ac]"1  ht_z 
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The  inverses  in  these  equations  may  be  shown  to  exist  by  construction  of 
jT  and  C_. 

Thus,  the  algebraic  approach  is  characterized  by  the  solution  of  a 
large  set  of  linear  equations.  For  example  in  the  EMI  scanner  [9]  the  un¬ 
known  imape  is  reconstructed  on  an  80  X  80  grid  which  gives  a  6400  equiva¬ 
lent  vector  dimension.  Also  =  180  projection  angles  are  used  with  m^  = 

160  points  per  projection  angle.  Thus,  K  is  approximately  28,000.  The 

A 

previous  solutions  for  M  if  directly  solved  would  require  inverting  6400  X 
6400  matrices. 

The  large  amount  of  computations  involved  motivate  two  avenues  of 
continued  study.  The  first  involves  the  use  of  iterative  solutions  of  the  large 
set  of  equations.  The  second  involves  studying  the  structure  of  the  H  matrix 
so  that  efficient  computational  algorithms  may  be  developed. 

Certain  structures  in  the  H  matrix  are  obvious.  For  example,  suppose 
the  coefficient  values  are  selected  proportional  to  the  intersection  area  of  a 
given  ray  and  the  elemental  unknown  area  as  shown  in  figure  4.  Then  assum¬ 
ing  either  one  or  two  rays  intersect  each  elemental  area  there  are  only  18 
possible  intersection  areas  as  shown  in  figure  la.  Thus,  although  the  inter¬ 
section  areas  changes  for  each  projection  angle,  the  formulas  for  the  com¬ 
putation  are  identical. 
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5.4  Application  of  Kalman  Filtering  to  the  Reconstruction  of  Images  from 
Noisy  Projections 
Firouz  Naderi,  Ali  Habibi 


A  digital  monochromatic  image  can  be  modeled  by  a  two  dimensional 
array  x(i,  j)  i ,  j  =  1,...,N  where  x(i,  j)  denotes  the  magnitude  of  the  gray 
level  at  point  (i,  j).  Then  a  one  dimensional  projection  of  the  NXN  image 
in  direction  0  is  defined  as  an  N  vector  YQ  where  each  component  of  Y_0  is 
the  sum  of  the  elements  of  the  image  along  a  line  making  an  angle  0  with  a 
fixed  reference  axis.  Figure  1  shows  examples  of  two  such  projections  for 
0  =  90  and  45  degrees.  Note  that  for  the  0  =  90°  projection 


where 


Y 

—90° 


[y 


l 

90°’ 


2 

y90°' 


yN  ]T 

y90°J 


Now,  let  x  be  an  N  vector  defined  as 


(1) 

(2) 


*  =  [x(l,  1 ),  x(l ,  2), 


x(l,N),x(2,  1),  . .  .x(N,N)]T 


(3) 


A  one  dimensional  projection  of  an  image  can  now  be  written  as 


—9  — 


(4) 
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2 

Where  the  N  by  N  matrix  CQ  is  called  the  projection  matrix  at  angle  6. 
For  the  example  of  figure  1  when  0=  90°  matrix  C.  is 


■^90° 


1  1  ...  1 

0  0  ...  0 

o  0  ...  0 


0  0  ...  0 

1  1  ...  1 

•  •  •  0 


.  0 

0  0...0  • . .  o 

1  1  . . .  1  0  ...  o 


0  1  1  ...  1 


Given  a  finite  number  of  projections  Y0  . Yg  and  the  corresponding 

projection  matrices  CQ  ,  . .  . . ^  it  is  cfesired  ^reconstruct  the  image, 
l.  e.  to  determine  the  N2  elements  x(i,  j),  i,  j  =  1 ,  . . .  ,  N.  This  problem 
can  be  viewed  as  an  attempt  to  solve  a  set  of  N  X  L  algebraic  equations  with 
N  unknowns.  In  most  cases  the  number  of  linearly  independent  equations 
available  is  not  equal  to  N  ;  therefore  obtaining  a  unique  solution  requires 
forming  a  meaningful  criterion  and  solving  the  equations  based  on  this 
additional  constraint. 


Another  approach  discussed  below  is  to  model  the  image  as  a  two 
dimensional  random  process,  and  use  some  a  priori  statistical  knowledge 
of  the  process  to  obtain  an  estimate  of  the  image. 

Reconstruction  of  Images  Using  Kalman  Filtering  The  approach  here 
is  to  derive  a  linear  dynamic  model  to  represent  the  imaging  system,  and 
then  use  linear  estimation  theory  to  obtain  an  estimate  of  the  image  which 
is  optimal  using  mean  square  error  criterion.  Let  Y (1 ),...,  Y(L)  be  L  one 
dimensional  projections  of  the  image  at  different  angles.  As  a  result  of 
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faulty  measurement,  Y(l),  . . . ,  Y(L)  are  not  the  true  projections,  but  contain 
some  noise.  Let  this  noise  be  an  additive  white  noise  denoted  by  the  N 
element  vector  n(k)  where  n(k)  is  the  noise  corresponding  to  projection 
_Y(k).  Then 

Y(k)  =  C(k)x  +  n(k)  (6) 

where  E[n(k)'^n(j)']  =  L6(k-j). 

Clearly  since  all  the  projections  correspond  to  the  same  fixed  image, 
the  dynamic  equations  representing  the  system  are 

x(k+l)  =  x(k)  (7a) 

Y(k)  =  C(k)x(k)  +  n(k)  (7b) 

Given  projections  Y_(l),  . .  .  ,  ^f(L)  it  is  desired  to  find  a  linear  estimate  of  x 
such  that  the  mean  square  error  is  minimum. 

The  recursive  solution  to  the  above  is  the  well  known  Kalman  estimator. 
Denoting  by  x(k+l)  the  estimate  of  the  image  obtained  after  observing  the 
projections  Y^(l), . . . ,  jf(k)  one  obtains 

x(k+l)  =  [l_  -  F(k)C(k)]x(k)  +  F(k)Y(k)  (8) 

where 

F(k)  =  P(k)TC(k)[C(k)P(k)TC(k)  +  LlJ'1  (9) 

and 

P(k)  =  E{[x(k)  -  x(k)][x(k)  -  x(k)]T}  (10) 

The  term  P(k)  is  the  covariance  of  the  error  after  k  iterations  and  has  the 
recursive  form 

P(k+1)  =  [l_  -  F(k)C(k)lP(k)  (11) 

Equations  (8),  (9)  and  (11)  provide  r  recursive  method  of  finding  the  estimate 
of  the  image.  However  along  with  this  equation  initial  conditions  for  x(k) 
and  for  I^(k)  are  needed. 
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A  reasonable  value  for  the  ini.'ial  conditions  are  [2] 

x(0)  =  E(x) 

P(0)  =  E[x  -  E(x)"l  [x  -  E(x)]T 

Some  aspects  of  this  approach  are  intriguing  enough  to  warrant  more 
investigation  of  this  method.  For  example: 

1)  The  error  covariance  equation  (eq.  11)  is  asymptotically  stable. 
Given  that  projection  matrices,  Cjk),  are  known  for  various  angles,  the 
error  covariance  can  be  calculated  for  all  k  before  any  projection  data  is 
actually  available.  This  suggests  that  it  might  be  possible  to  decide  how 
many  projections  are  necessary  to  achieve  a  certain  fidelity. 

2)  The  error  covariance  depends  on  the  projection  matrices  which 

in  turn  depend  on  the  angle  of  projection.  It  might  be  possible  then  to  decide 
the  angles  of  projections  more  judiciously  by  proper  evaluation  of  the  error 
covariance. 

At  the  present  the  biggest  disadvantage  of  the  method  is  the  excessive 
computations  which  are  required  because  of  using  large  matrices.  However 
this  method  has  distinct  advantages  over  other  methods  that  are  used  for 
reconstructing  images  from  noisy  projections.  Namely,  it  is  iterative  and 
provides  a  measure  for  the  error  and  is  based  on  an  established  estimation 
theory  procedure. 
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5.  5  Boundary  Estimation  of  Statistical  Objects  in  Noise 
Nasser  E.  Nahi  and  Mohammad  Jahanshahi 

The  problem  considered  is  determination  of  the  boundary  of  an  object 
in  noise.  Standard  boundary  detection  techniques,  in  general,  utilize  some 
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form  of  differencing  or  differentiation,  [l]  .  These  methods,  however,  are 
sensitive  to  noise.  Various  extensions  of  these  techniques  appear  in  [2], 
[3]  .  They  are  mainly  effective  for  high  signal  to  noise  ratio  images. 

The  images  considered  in  this  work  can  be  partitioned  into  two 
regions:  background  and  foreground.  The  foreground  is  assumed  to  form 
a  horizontally  convex  set  in  the  plane. 

Definition:  A  set  E  CR  is  said  to  be  horizontally  convex  if  given 

—  =  (x1»x2^eE*  £  =  (Yj*  V2)*:E  with  x ^  i  y^  and  x ^  =  y 2>  then 

wx  +  (1  -w)^eE,  where  0  <  w  <  1. 

Example:  Sets  E^  and  E^,  in  figure  1,  are  horizontally  convex. 

Set  E^  is  not. 

An  object  of  interest  whose  intensity  levels  dominate  those  of  the 
background  will  always  be  assumed  to  exist  in  the  image.  This  object 
forms  the  foreground. 

Problem  Statement  Assuming 

1*  image  consisting  of  a  background  and  a  horizontally  convex 

object  is  given; 

2.  The  image  is  represented  by  b(z,  n)T(z,  n),  where  b(z,  n)  is  a 

random  process  representing  grey  levels  within  the  object,  and 

!1  when  (z,  n)  is  a  point  in  the  object, 

0  when  (z,  n)  is  a  point  in  the  background; 

3.  The  process  b(z,  n)  is  characterized  in  a  statistical  sense,  in 
terms  of  its  first  ana  second  order  moments; 

4.  The  operator  T(z,  n)  def  ned  the  object  boundary; 

5.  The  observables 

y(z,n)  =  b(z,  n)T(z,  n)  +  v(z,  n)  , 

where  v(z,  n)  represents  the  observation  noise  with  known  first 
and  second  order  statistics,  are  given. 


The  goal  is  to  find  an  estimate  of  the  object  boundary. 


A  Solution  Method  The  image,  as  stated,  is  characterized  by  a  two 


dimensional  random  process.  Since  the  solution  of  the  proposed  method 
will  utilize  conventional  digital  computers  with  litter  or  no  parallel-computa¬ 
tion  capabilities,  an  image  scanner  is  employed.  Therefore,  the  problem  is 
restated  as  the  following:  Assuming 

1.  The  observation 


y(t)  =  s  (t)X(t)  +  v(t) 


for  0  st  s  T1 


is  given,  where  s(t)  denotes  the  output  of  a  line  scanner,  represent¬ 
ing  the  grey  level  of  the  image  at  time  t,  and  where 


m. 


with 


X(t)  =  u(t-a^)  -  u(t-c^) 

A-m  ^ 


m1  =  first  line  of  the  image  containing  the  object; 
m2  =  last  line  of  the  image  containing  the  object; 
a.  =  start  of  the  object  in  line  A ; 

Cg  =  end  of  the  object  in  line  A  ; 
u( • )  =  step  function; 


and 


v(t)  =  process  representing  the  inaccuracies 
introduced  by  the  scanner  or  any  other 
phenomenon  distorting  the  image. 

2.  The  sequence  x  =  (a^  ,  c^  ),  A  =  m^,  m^  +  1,  .  . . ,  m^,  is  Markov. 

3.  The  density  function  f(Xg  |xg  is  Gaussian. 

Find  the  M.A.P.  estimates  of  a^  and  c^  ,  m^  ^  l  im2'  denoted  by  a^  and 
,  respectively,  and  determine  the  values  of  m^  and 

A  recursive  solution  to  the  above  problem  has  been  obtained.  Figure 
2  contains  two  examples  showing  the  feasibility  of  the  boundary  estimation 
technique.  Further  investigations  are  in  progreLS  as  to  the  refinement  of 
some  of  the  theoretical  aspects  of  the  problem. 


(a)  original  diamond 


(c)  noisy  diamond 


(d)  noisy  square 


,(e)  boundary  estimate 


(£)  boundary  estimate 


(b)  original  square 


Figure  5.  5-2.  Boundary  estimation  examples. 
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6.  Image  Analysis  Projects 

The  image  analysis  projects  are  concerned  with  the  background 
technology  necessary  to  effectively  design  image  coding,  restortion, 
enhancement,  and  data  extraction  systems.  Of  particular  interest  are 
models  of  the  human  visual  system  for  monochrome  and  color  images, 
and  the  development  of  quantitative  measures  of  image  fidelity  and 
intelligibility. 

The  following  project  traces  the  development  of  a  model  of  human 
color  vision.  This  model  has  been  found  to  accurately  predict  known 
visual  phenomena  such  as  spatial  color  contrast  and  color  blindness. 

6.  1  A  New  Modei  of  Color  Vision  and  Some  Practical  Implications 

Werner  Frei 

Color  image  coding  and  processing  techniques  usually  operate  on 
three  scalar  functions  of  some  space  variables  such  as  red,  green  and 
blue  color  separations.  Under  the  laws  of  colorimetry,  these  scalars 
called  tristimulus  values  represent  light  energy  quantities  and  the  operations 
defined  on  them  by  color  matching  experiments  conveniently  satisfy  the 
requirements  of  linear  mappings;  for  example,  the  additive  mixture  of 
two  arbitrary  colors  is  mapped  into  the  sum  of  their  respective  tristimulus 
value  s . 

In  terms  of  the  human  observer,  however,  colorimetric  color  repre¬ 
sentations  have  serious  disadvantages,  because  the  tristimuli  are  not 
representative  of  the  physiological  response  evoqued  in  the  human  visual 
system.  It  has  been  shown  that  the  visual  response  is  a  highly  nonlinear 
function  of  the  energetic  tristimulus  quantities.  Tristimulus  spaces  are 
therefore  not  euclidean,  and  unfortunately,  ordinary  superposition  of 
tristimuli  does  not  entail  superposition  of  responses.  For  example,  if 
the  tristimulus  values  of  a  color  are  doubled,  one  does  not  obtain  a  color 
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that  appears  twice  as  bright,  twice  as  saturated,  etc.  This  remark  applies 
of  course  to  spatial  arrangements  of  colored  stimuli  so  that  visual  modula¬ 
tion  transfer  functions  are,  strictly  speaking,  not  defined  in  tristimulus 
space. 

It  is  therefore  no  surprise  that  distortion  measures  such  as  mean 
square  error,  or  visual  modulation  transfer  functions,  etc.  are  of  little 
relevance  to  the  observer  of  a  picture  processed  under  that  kind  of  criterion. 

In  view  of  these  difficulties,  the  simple  visual  model  developed  here 
has  a  variety  of  practical  implications  in  image  processing  and  possibly  in 
vision  research.  It  is  shown  that  its  structure  enables  one  to  define  an  alge¬ 
braic  system  in  tristimulus  space  with  a  generalized  superposition  consis¬ 
tent  with  known  perceptual  phenomena. 

As  a  result,  one  can  define  error  measures  in  a  euclidian  domain 
that  is  a  linear  vector  mapping  of  tristimulus  space,  given  a  particular 
set  of  basis  vectors.  This  gives  analytical  convenience  in  the  design  and 
evaluation  of  coding  and  processing  techniques. 

On  the  other  hand,  with  superposition  defined,  visual  MTF  can  be 
defined  ard  it  can  be  shown  that  spatial  color  contrast  phenomena  such  as 
sim  iltaneous  color  contrast,  color  shadows,  color  mach  bands,  etc.,  can 
be  modeled  by  linear  spatial  filters. 

Superposition  is  also  expected  to  enhance  the  effectiveness  of  pseudo¬ 
color  techniques.  Since  the  visual  system  is  able  to  recognize  more  or 
less  independently  three  attributes  of  colors  (such  as  brightness,  hue  and 
saturation),  pseudo-color  techniques  attempt  to  display  multiple  functions 
(at  most  three)  of  the  space  variables  as  a  single  composite  color  picture. 

If  the  observer  is  expected  to  be  able  to  recognize  the  indi\ idual  functions 
in  the  pseudo-color  image,  it  is  quite  evident  that  the  functions  of  interest 
must  be  mapped  into  color  space  in  accordance  with  subjective  superposition. 


The  visual  model  discussed  here  provides  a  tool  to  design  precisely  that 
kind  of  mapping,  and,  in  addition,  it  predicts  the  spatial  resolution  to  be 
expected  for  each  function  mapped  into  color  space. 


Visual  Model  Figure  1  shows  a  block  diag  ram  of  the  model.  The 
basic  ideas  and  the  physiological  evidence  justifying  its  structure  have 
been  discussed  in  a  previous  report  [l].  In  brief,  the  first  stage  converts 
the  input  spectral  energy  distirbution  of  light  C(X)  (a  function  of  the  space 
variables  x  and  y)  into  a  tristimulus  vector  t  =  t(x,y)  whose  three  components 
represent  the  respective  amounts  of  light  energy  absorbed  by  three  types  of 
photo-receptors  (Young-Helmholm  theory): 

t  =  Cvt^f  tj  =  t.(x,  y) 

XTT 

t.  =  I  C(X)t.(X)dX  i  =  1,2,3 


where  t.(X)  is  the  spectral  absorbtion  function  of  the  i-th  type  of  receptor. 
The  second  stage  of  the  model  represents  the  conversion  of  absorbed 
energies  to  neural  signals,  according  to  an  approximate  logarithmic 
relationship  (Weber-Fechner  type  of  response).  The  components  of  the 
resulting  vector  t*  are:  log  t^,  log  log  t^  (the  superscript  star  denotes 
the  logarithmic  domain).  The  last  stage  of  the  model  contains  a  matrix  of 
spatial  filters  with  transfer  function 


H*(ju*,jv*)  = 


r  H* 


■  h; 


l-h? 
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0 
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2 
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H 


0 
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3  J 


Hii  = 

J 


where  u*,  v*  are  spatial  frequencies  in  the  logarithmic  domain.  This 
stage  simultaneously  represents  two  phenomena; 

(a)  linear  differences  (inhibitions)  between  the  outputs  of  the  two  types  of 
receptors  (demonstrated  in  the  retinas  of  monkeys  by  deValois  [2]).  These 
differences  are  represented  by  equal  terms  of  opposite  sign  within  a  row  of 
H*.  They  can  be  seen  to  generate  two  "  perceptual"  chromatic  signals  by 
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Figure  6.  1-1.  Model  of  color  vision 


observing  that 


4  -  H2(t2  -  tt) 


g*  =  H*(t*  -  t*) 


& 


Under  the  assumption  of  exact  logarithmic  receptor  responses,  g*  and  g* 

U  J 

are  functions  of  the  ratios  of  pairs  of  tristimulus  values,  e.  g.  chromatic 
quantities  independent  of  the  absolute  light  intensity.  (This  approximation 
is  not  valid  at  very  low  or  very  high  intensities,  but  quite  acceptable  within 
most  of  the  range  of  interest.  The  intensity  is  represented  by 

J/ 

g*  =  logft^ni 

(b)  The  matrix  H*  also  represents  weighted  linear  summations  between 
the  outputs  of  spatially  distributed  like  receptors.  Assuming  shift -invariance, 
suppose  that  each  summation  output  g*(x.,y.)  is  related  to  the  outputs  of  the 
receptor  array  t*(x^,  y^  )  by  the  sum 

=  ^  £t*<VVh*(Wyj-V 

where  h*(x,  y)  denotes  the  impulse  response  matrix  in  the  logarithmic 
domain.  Taking  the  discrete  two-dimensional  Fourier  transform,  gives 


g.(x,  y)  =  7 


^  jC  y>]H^(ju*.  jv*^ 


(Note:  the  shorthand  notation  gT  =  Y  H]'.t'.'  is  used  instead  of  the  above 

1  j  U  J 

expression. )  Letting  t^(X)  be  equal  to  the  luminous  efficiency  function 
V^,  and  considering  the  scalar  expression  for  ,  a  familiar  model  of 
achromatic  vision 


,'[(x,y)  =  7  ^[log  I(x,  y)’|H"'(ju*,  jv*)) 
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is  obtained  [3*].  This  model  has  been  successfully  used  to  predict  various 
visual  phenomena  such  as  Mach  bands,  simultaneous  contrast,  etc.,  as 
well  as  for  image  coding  and  processing  [4]. 

Spatial  Color  Contrast  Phenomena  Examination  will  now  be  directed 
toward  the  mechanisms  by  which  chromatic  contrast  effects  such  as  color 
shadows,  etc.  ,  can  be  represented  by  the  model.  Suppose  that  a  small  test 
spot  with  tristimulus  values  tj,t  ,t3  is  viewed  against  a  colored  background 
with  values  *10*  *20' *30*  Assuming  total  lateral  inhibition,  e.  g.  mutual 
subtractions  of  the  outputs  from  neighboring  receptors,  the  apparent 
chromatic  quantities  g*',  g*'  and  g*^,  of  the  test  color  and  background 
respectively,  become 


In  other  words,  the  background  appears  to  be  white,  while  the  test  color 
is  shifted  towards  a  color  complementary  to  that  of  the  background.  The 
same  apparent  chromaticities  are  predicted  by  the  model  if  one  assumes 
that  the  responses  of  the  photo- receptors  are  decreased  in  proportion  to 
the  average  energies  absorbed  (bleaching).  The  new  tristimuli  are 

t. 

t!  =  -L 
1  ti0 

And  the  apparent  chromatic,  quantities  become 
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The  apparent  chromatic  shifts  predicted  correspond  to  the  vonKries 
coefficient  law  regardless  of  the  physiological  assumption  we  make 
(inhibition  or  bleaching). 

In  practice  however,  the  visual  system  does  not  ignore  completely 
a  non-white  average  sceme  chromaticity,  particularly  at  saturated  average 
chromaticities.  Letting  lateral  inhibition  be  weighted  by  a  factor  k*  <  1 
the  following  apparent  chromaticities  are  obtained 


•  lr)  -  •;  '*■  (f) 


k*  * 1 


k2  <  1 


Since  the  influence  of  various  colored  areas  within  the  visual 
field  on  each  other's  corresponding  perceptual  quantities  is  a  function 
of  their  geometric  distance,  let  k|,  k^'  be  functions  of  the  geometric 
d’ stance  between  the  excited  receptors  and  write  for  each  summation  output 

t2(Xk’y/) 

g2(VV  *  7  Vvwv 

51!  ^3^Xk,y/^ 

g3(x  ,y  )  =  22/2— - —  k2(Xi"Vyi"y/) 

3  k  l  tj(xk,y/)  2  1  *  1  1 

It  is  then  observed  that  k  (x,y)  is  an  impulse  response,  and  the  color 
contrast  phenomena  can  be  represented  by  the  spatial  filters  H*(ju*,  jv*), 
H*(ju*,jv*) 
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Figure  2  shows  a  comparison  of  measured  visual  modulation  transfer  functions. 
The  question  whether  the  chromatic  response  is  attenuated  at  very  low  spatial 
frequencies  (spatial  band-pass  response)  is  subject  to  controversy.  Note 
that  the  measurements  of  chromatic  responses  were  not  made  in  a  logarithmic 
domain. 

Spectral  Sensitivities  of  the  Receptors  According  to  Konig's  theory, 
color  blindness  known  as  dichromatism  occurs  whenever  one  of  the  three 
types  of  photo-receptors  is  inoperative  in  some  fashion.  The  theory  there¬ 
fore  predicts  the  existence  of  the  three  types  of  dichromats  whose  color 
matches  are  uniquely  determined  by  two  of  the  three  energy  absorption 
quantities 


t.  =  J  S(X)t.(X)d(X)  i  =  1,2 


Under  Konig's  assumption,  colors  with  two  fixed  values  t.,L  and  arbitrary 
value  t^  should  all  look  alike  to  the  corresponding  dichromat.  It  can  be 
shown  that  the  chromaticity  of  such  equivalent  colors  are  sets  of  converging 
straight  lines  in  the  CIE  x-y  chromaticity  diagram,  and  that  the  centers  of 
convergence  are  the  chromaticities  of  the  respective  missing  primaries. 


Experiments  show  that  there  exist  indeed  three  types  of  dichromats, 
called  protanopes,  deuteranopes  and  tritanopes  and  that  their  color  matches 
correspond  well  with  the  predicted  behavior  (fig.  3). 


Let  A  be  a  3  x  3  matrix  which  maps  the  CIE  XYZ  tristimulus  space 
into  the  model's  tristimulus  space 
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CHROMATIC  RESPONSE 
(von  der  Horst) 


PROTANOPE  DEUTERANOPE  TRITANOPE 
loci  of  dichromatic  color  matches 


Figure  6. 1-3.  Confuaion  loci  for  the  three  types  of  dichromats  in  the 
CIE  x-y  chromaticity  diagram. 
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Also  let  x^,  y^>  z^»  k  =  1,  ...  3  be  the  chromaticity  coordinates  of  the 
convergence  for  the  protanope,  deuteranope  and  tritanope  respectively. 
Since  x^.y^  zy  are  chromaticities  of  the  unknown  primarties,  the 
following  equations  can  be  written  to  determine  A  =  [a..] 


Three  more  equations  can  then  be  written  by  defining  an  arbitrary  reference 
white  for  the  model,  for  which  tj=t  =t^=l.  Finally,  the  spectral  sensitivity 
functions  of  the  receptors  are  given  by 


‘i<M  *x 

t2(X)  =  [at  7x 

t3(x)J  Ux. 


where  x  ,  y  ,  z  are  the  CIE  color  matching  curves  for  the  standard  observer. 

A,  A  A, 

Using  essentially  the  method  outlined,  and  rounding  Pitts  coordinates  for  the 
deuteranopic  convergence  point  from  x  =  1.08,  y  =  -0.08  to  x  =  1.00  and 
y  =  0.  0,  Judd  obtained  the  following  transformation  matrix 


A  =  -0.460  1.359  0.  101 


The  spectral  sensitivities  of  the  receptors  thus  obtained  are  unimodal 
positive  functions  of  the  wavelength,  as  one  would  except  from  a  physiological 
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standpoint  (fig.  4).  Furthermore,  one  of  the  receptors  exhibits  a  V  ,  or 

X 

luminance  response,  as  is  required  for  our  model  to  be  considered  with 
Abney's  law  of  luminance  addition.  Figure  4  also  shows  the  remarkable 
agreement  between  the  measured  luminous  efficiency  functions  of  the  normal 
and  dichromatic  observers,  and  the  luminous  efficiencies  predicted  with 
the  model. 


Uniform  Perception  Space  g*  While  the  visual  model  has  been  based 
upon  physiological  evidence,  quite  a  few  approximations  have  been  made.  It 
remains  to  verify  whether  the  model  predicts  measured  color  differences. 

First  note  that  g*  =  log  luminance  is  in  agreement  with  Weber's  law 

of  perception  for  intensities.  As  far  as  the  chromatic  quantities  g*  and  g* 

2  3 

are  concerned,  consider  concentric  circles  centered  at  the  origin  and 
straight  lines  radiating  from  the  origin  in  the  g*,  g*  plane.  If  the  model 
were  perfect,  the  circles  would  be  loci  of  constant  saturation  colors,  and 
the  lines  loci  of  constant  hue  colors.  Figure  5  shows  a  set  of  such  equi¬ 
distant  circles  and  lines,  mapped  into  the  CIE  x-y  chromaticity  diagram. 

The  agreement  with  the  Munsell  equi-distant  set  is  remarkable  if  scaling 
constants  Cj,c2  (to  be  included  in  the  filter  function  H*  and  H*)  are  introduced 
such  that 

*  .  *2 
e2  =  cj  log- 


„*  - 


Cl  lQg  ~ 


"l 


Assuming  that  the  sensation  of  brightness  is  evoked  by  the  total  activity  of 
the  three  perceptual  channels,  brightness  is  then  equal  to  the  norm  of  g*, 

|g*l  =  (gj  +  g'2  +  g^  )*•  The  mapped  set  of  circles  and  lines  now  repre¬ 

sents  equi-brightness  geodesics.  A  remarkable  agreement  with  the  geodesics 
obtained  with  Stile's  line  element  is  noted  [5"],  Improvements  may  be  expected 
by  replacing  the  logarithmic  approximation  by  a  more  realistic  function  of 
the  type  AI/K+I.  However,  the  original  approximation  shall  be  retained  for  the 
next  discussion. 


CIE  V^-curve 
Konig  g,  (Judd) 
Vx  proto  nope 


Comparison  of  Konig-type  receptor  sensitivities  and  luminou 
efficiency  functions  of  normal  and  dichromatic  observers. 


CIE  Spectrum  locus 


Figure  6,1-5,  Equi-distant  loci  of  constant  hue  and  constant 
saturation  Dredicted  bv  the  model. 


Generalized  Superposition  in  Tristimulus  Space  Define  two  vector 


spaces  V  and  W  over  the  field  of  real  numbers  reR,  with  teV,  t  =  (tj,t  ,t  ) 
and  g  eW,  g*  =  (gj\  g|,  gp.  Let  the  components  of  t  be  tristimulus 

values  with  respect  to  the  spectral  receptor  sensitivities  of  the  model  and 
some  arbitrary  reference  white,  and  define  the  following  laws  of  composition 
on  V 


Wv 


*A  +  lB  "  (tlA*lB't2At2B’t3At3B) 

r  r  r 

rt  =  (tj,t2,t3) 


Next,  let  the  components  gj,  g£,  of  g*  be  the  perceptual  quantities 

defined  in  the  model,  and  let  ordinary  addition  and  multiplication  by  scalars 

be  the  laws  of  composition  on  W.  Since  the  three  receptor  sensitivity  curves 

selected  are  positive  between  X  ,  X  (the  limits  of  the  visible  spectrum), 

J-/  u 

*1**2' *3  arC  always  larger  or  equal  to  zero.  Excluding  zero  light  intensities 
(which  is  quite  acceptable  for  practical  images)  the  vector  mapping 


T:  t 


g 


* 


defined  by  the  model  as 


^1'  *2'  *3^ 


/  f  -  l2  *3  \ 

(H1  1o8*1-H2  1o*^-H3  l°*rJ 


can  easily  be  shown  to  be  linear.  The  inverse  mapping  T  is 

T  1t  (gj,  g'2,  gp - ♦  |exp(gp1/Hi,  exp(g*,  g|)1/H^,  exp(g*,  g*)1/H3^ 


It  is  seen  that  the  operations  defined  on  t  (which  is  an  element  of  tristimulus 
space,  relative  to  a  particular  basis  determined  by  the  model's  receptor 
sensitivities)  correspond  to  ordinary  addition  and  multiplication  by  scalars 
of  perceptual  quantities  in  the  euclidean  space  W.  For  example,  multiplying 
the  tristimuli  of  a  color  by  a  scalar 
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r(trW  =  *  (rgi’ rg2' rg3) 

results  in  a  color  which  appears  r-times  as  light,  r-times  as  saturated  and, 
letting  brightness  be  approximately  proportional  to  the  norm  of  g*,  a  color 
r-times  as  bright  as  the  original  color  is  obtained. 

Conclusions  Generalized  superposition  thus  defined  opens  a  broad 
variety  of  image  coding  and  processing  related  problems  to  analytical 
treatment.  In  particular,  homomorphic  filtering  can  be  readily  applied 
to  color  images,  if  tristimuli  are  referred  to  the  basis  indicated.  The 
simple  mapping  T:  t  — h >g*  can  be  used  in  any  real-time  environment  to 
solve  the  old  prob]em  of  optimal  color  quantization.  Finally,  it  is  hoped 
that  the  reflections  on  superposition  will  help  to  clarify  some  of  the 
controversy  on  chromatic  modulation  transfer  functions  in  the  visual 
system. 
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7. 


Image  Processing  Hardware  and  Software  Projects 


The  image  processing  hardware  and  software  projects  are  devel¬ 
opmental  projects  supportive  of  the  image  processing  research. 

The  first  project  describes  the  progress  toward  the  development 
of  a  real  time  digital  image  display  system  for  the  ARPANET.  This 
device  which  is  connected  to  a  port  on  the  TIP  provides  a  flicker  free 
display  of  digital  images  transmitted  over  the  network. 

In  the  second  report  plans  are  outlined  for  the  construction  of  a 
digital  image  recording  and  display  system.  This  system  will  be  capable 
of  recording  real  time  color  television  for  playback  into  a  digital  computer; 
and  will  record  computer  generated  images  for  real  time  playback. 

The  last  report  is  concerned  with  progress  on  software  development 
for  image  processing  and  networking. 


7.1  Real  Time  ARPANET  Digital  Image  Display 
Toyone  Mayeda 

The  digital  color  image  display  for  use  on  the  ARPANET  is  presently 
ready  for  software  checkout.  The  refresh  memory  integrated  circuits  were 
finally  received  during  February  1974. 

A  block  diagram  of  the  system  is  shown  in  Figure  1,  and  photographs 
of  the  system  are  shown  in  Figure  2.  The  display  system  specifications 
are  listed  below: 

1.  Receive  asynchronous  digital  picture  information  from  the 
ARPANET  TIP  with  brightness  resolution  up  to  64  levels  and 
at  input  rates  up  to  17.  2K  band; 

2.  Store  the  received  data  in  an  array  of  up  to  256  X  256  six  bit 
picture  values; 

3.  Display  a  true  6  bit  monochrome  image  on  the  monitor; 

4.  Display  a  pseudo  color  image  by  use  of  a  random  access  memory 
which  is  addressed  by  *;he  refresh  memory  output.  The  RAM  can 
be  remotely  programmed  from  the  TIP  or  by  local  switch  con- 


INPUT 


igure  7. 1-1.  ARPANET  digital  image  display. 


(b)  total  system  -  RCA  color  TV,  processor  electronics 
and  power  supply. 


Figure  7.  1-2.  ARPANET  digital  image  display 
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trol.  Over  4096  different  color  combinations  of  hue,  saturation 
and  luminance  are  available  for  pseudo-coloring. 

Design  of  an  8  bit  monochrome  digital  image  display  has  been 
completed  and  is  presently  in  the  wire  listing  phase  (prior  to  fabrication). 
The  specifications  of  this  system  are: 

1.  Receive  asynchronous  digital  picture  information  from  the 
ARPANET  TIP  with  brightness  resolution  up  to  256  levels 
and  at  input  rates  up  to  19.  2K  baud; 

Include  a  function  memory  which  can  be  used  to  translate  the 
8  bit  data  words  (from  the  refresh  memory)  with  any  desired 
transfer  curve.  The  function  memory  can  be  remotely  pro¬ 
grammed  from  the  TIP  or  by  local  switch  control; 

Display  a  256  x  256  eight  bit  image,  or  a  six  or  seven  bit  image 
and  a  one  bit  graphic  overlay. 

Use  an  alphanumeric  keyboard  to  communicate  with  the 
ARPANET  T TP  and  also  to  generate  alphanumeric  char¬ 
acters  on  the  display  monitor. 

Output  the  monochrome  video  data  and  alphanumeric  charac¬ 
ters  in  composite  RF  format  so  that  it  can  be  displayed  on  any 
TV  receiver  using  its  antenna  input. 


3. 


4. 


5. 


7.  2  Real  Time  Color  Image  Digital  Recorder  and  Display 
Toyone  Mayeda 

Design  has  been  started  on  a  real  time  image  digital  recorder  and 
display  system  as  shown  in  Figure  1.  The  system  will  initially  be  designed 
to  process  an  eight  bit  monochrome  image,  but  eventually,  will  have  the 
capability  to  process  a  color  image  using  additional  tracks  of  the  recorder. 
The  digital  recorder  and  analog  to  digital  converter  have  been  ordered. 

With  reference  to  the  digital  recorder  in  Figure  1,  the  monochrome  system 
operates  in  the  following  four  modes: 

1.  In  the  600  ips  record  mode,  video  from  the  monochrome  TV 

camera  is  digitized  at  approximately  10  MHz  by  the  A/D  conver- 
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Figure  7*2-1.  Realtime  image  digital  recorder  and  display. 
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ter  and  the  eight  digitized  bits  are  selected.  A  ninth 
channel,  which  consists  of  a  continuous  16  bit  sync 
word  pattern,  is  distributed  with  the  eight  data  channels 
and  all  nine  channels  are  recorded. 

2.  In  the  17/8  ips  playback  mode,  the  nine  channels  are  aligned 
in  the  deskew  buffer.  After  the  sync  words  are  separated,  the 
data  in  the  eight  remaining  channels  are  transferred  to  the 
controller  to  be  processed. 

3.  In  the  17/8  ips  record  mode,  the  processed  data  from  the 
controller  is  recorded  with  the  sync  word  as  in  mode  1. 

4.  In  the  600  ips  playback  mode,  the  nine  channels  are  aligned 
and  the  sync  words  are  separated  as  in  mode  2.  The  eight 
remaining  data  channels  are  applied  to  a  D/A  converter  and 
the  output  video  data  is  displayed  on  a  TV  monitor. 

7.  3  Software  Progress  Report 
James  Pepin 

In  the  last  six  months  several  important  projects  have  been  initiated 
by  the  programming  staff.  The  most  visible  is  a  new  network  timesharing 
monitor.  The  second  area  of  work  is  the  interfacing  of  the  Optronics 
scanner  and  digitizer.  A  third  area  of  effort  has  been  in  the  front-end 
image  processing  system  design.  Finally  a  PDP-11/10  has  been  added  to 
the  Image  Processing  Instite  (IPI)  hardware.  This  machine  well  act  as  a 
front-end  to  the  HP2100,  relieving  it  of  some  device  control  functions. 

The  new  netv/ork  support  has  many  invisible  improvements  as  well 
as  some  that  are  at  once  apparant.  A  password  capability  has  been  added 
to  allow  users  data  integrity.  This  also  will  allow  more  awareness  of  who 
is  using  the  system  and  what  they  are  doing.  Another  important  feature 
added  to  the  monitor  is  the  capability  to  perform  multi-tasking  within  the 
timesharing  system  itself.  This  will  allow  implementation  of  some  compli¬ 
cated  network  protocals  without  complex  code  in  each  program.  The  system 
also  has  been  made  more  efficient  in  its  input/ output  and  user  swapping 
areas.  This  improvement  should  be  very  noticeable  in  heavy  loading  situa- 
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tions.  Other  improvements  in  this  area  include  new  and  better  help  files. 

A  file  structure  scheme  has  been  implemented  that  will  allow  for  faster  and 
easier  searches  of  HELP  and  other  similar  data  bases. 

In  December  IPI  received  an  Optronics  microdensitomiter.  This 
device  has  required  a  large  software  effort  for  interfacing  to  the  HP2100. 
This  microdensitometer  has  maiy  options  that  must  be  implemented  by 
the  computer  so  that  a  user  can  scan  or  display  an  image  without  detailed 
knowledge  of  the  control  logic. 

At  the  end  of  January  the  Institute  received  a  PDP  11/10.  This 
computer  will  be  used  as  a  device  controller  to  relieve  the  HP2100  of  some 
of  its  controlling  duties.  As  it  now  stands  the  2100  has  to  handle  all  the 
'bit  watching'  for  the  devices  in  the  lab;  with  the  addition  of  the  11/10  some  of 
this  task  can  be  moved  to  the  11/10,  allowing  the  2100  to  be  free  to  handle 

other  tasks.  This  technique  will  enable  operation  of  two  or  more  devices 
simultaneously. 

During  the  last  six  months  the  staff  has  been  participating  in  several 
committees  that  have  as  their  responsibility  the  formulation  of  the  needs  of 
the  image  processing  community  on  the  ARPA  net.  There  have  been  two 
areas  of  prime  interest:  software  requirements  and  hardware  requirements. 
The  software  area  investigation  has  lead  to  the  conclusion  that  some  sort  of 
insulation  is  necessary  between  the  average  image  processing  user  and  the 
large  machines  on  the  net.  The  approach  taken  is  that  a  software  package 
should  be  implemented  that  will  present  to  the  image  processing  user  one 
command  language.  This  will  allow  the  user  to  use  any  large  machine 
while  only  learning  the  command  for  image  processing.  Following  this  line 
of  thought  it  follows  that  it  would  be  desirable  to  have  an  image  processing 
'front-end'.  This  would  be  a  computer  system  that  would  run  this  common 
language,  convert  that  into  the  Job  Control  of  the  proper  machine,  and  then 
send  the  JCL  along  with  all  pictures  and  other  data  needed  to  complete  the 
user’s  task.  This  system  could  also  act  as  a  conferencing  tool  for  image 
processors  and  could  be  a  common  collection  point  for  literature  in  the  field. 
Various  hardware  and  software  approaches  have  been  investigated. 


-150- 


Pub  1 1 ca  1 1 ons 


resulted  from  ARPA  sponsored  research. 


Harry  C.  Andrews,  "Digital  Fourier  Transforms 

!S?!"epP  m-ur"  App"art  °PtiCS'  VoK  13' 


as  a  Means  of 
No.  1,  January, 


Harry  C.  Andrews  and  C.  K.  Chen,  "Nonlinear  Intrinsic  Dimen- 
C-2?P  No!'  2rSh?ii?“'nnfE  TranSaCtions  °"  Computer,  Vo.. 

&eV5ors7;  ;rn^,r^??storat,pn-survay<" 

;s  ”*'■*  c!-“ 


L.  D.  Davisson,  "Universal  Noiseless  PnHimx  "  i  err  t--,.. 
on  Information  Theory,  November,  1973  R/  E  Transact,ons 

^«epthL*  HaH/  "Almost  Uniform  Distributions  for  Computer  Im¬ 
age  Enhancement,  IEEE  Transactions  on  Computer  s,  February,  1974 . 

Ernest  L.  Hall  and  William  0.  Crawford,  "Automated  Class!  f|  ca- 

26th  ^  r?  ?n  °PacItles  In  Radiographs  of  Coal  Workers  " 

September3  1973fefe?Ce  °n  EngineerinK  In  Medicine  and  Biology; 
eptember,  1973.  (also  supported  by  NIOSH  grant  HSM-99-72- 30.) 

Ernest  L.  Hall,  Kendell  Preston,  and  F.  E.  Roberts  "Clascffi- 

F I  rstn|nternat  lin10?  ?f  ?lack  Lung  Disease  from  Chest  X-rays," 
1Q7<  JnfRrnat,onal  Conference  on  Pattern  Recognition,  October 
1973.  (also  supported  by  NIOSH  grant  HSM-99-72- 30. ) 

E.  L.  Hall,  R.  p,  Kruger  and  A.  F.  Turner,  "An  Optical  niri^i 
PubH^^TpUc^  Engineer "ng?^  ^  X-rayS'"  f°r 

*  - '  an?„pop!  ;?r0f  f  i™- ir  sn£s 

wn.'Tt.^tarLrMr'p”!-?^?:5  c°nferencp- 


-151- 


All  Habibi,  "Hybrid  Coding  of  Pictorial  Data,"  IEEE  Transac¬ 
tions  on  Communications,  Vo  1 .  COM-22,  No.  5,  May,  19  74. 

D PCm' 'a nH b,Tr -»ROfa  1  d  "Rrshe1'  "A  Unified  Representation  of 

n  i  ca t?ons JrVol  f°COM-2  2,  "wof^^May" '  1974yranSaCt  *  °nS  “ 

rnd|tnf'i'b!crSTr  S*  R°b,nson'  "A  Survey  of  Digital  Picture 
Coding,  i EEE  Computer (special  i ssue  on  Dig! tal  ImageProcess- 

?!rl0H?mibi/DWiniam  K<  Pratt'  Guner  Icbinson,  et.al.,  "Real 
ninuesm"Eiq7Udyn^anCy  Ref,uft,°n  Usinp,  Transform  Coding  Tech- 
1974?  '  1974  ,nternatlona1  Communications  Conference,  June, 

D PCM  rni?'  ",l?a?e  f1°de11,nK  for  Unification  of  Transform  and 

Se?°edien?  Uatl°nal  E1*Ct™"S 

Anil  K.  Jain,  "Computational  Considerations  in  Generation  of 
Me°t?nLCS1^3  SpaCe'"  0Pt'Ci,,  S°ciety  of  ^erfca  Annuaf 

and 1  aJd  Edward  Angel,  "Image  Restoration,  Modelling 

and  Reduction  °f  Dimensionality,"  accepted  for  pub  1 i ca t i on  i n 
IEEE  Transactions  on  Computers,  May,  1974.  "cation  in 

"Comnuter ' i  Ml]1ian  R*  Tho"Pson,  and  A.  Franklin  Turner, 

Systems^an  and°Cybernet I*cs^1°JanuaryS#1974EE(a^soSaCt '°n.  S"k 

Ml  OSH  grant  HSM-99-72-30. )  '  1974*  (also  su™orted  by 

Richard  P.  Kruger,  Ernest  L.  Hall,  and  A.  Franklin  Turner  "a 
Screening  System  for  Pneumoconiosis  Chest  X-Ray 

ics  Bo^fAn  m  *  mE  Colnference  on  Systems  Man  and  Cybernet- 
graAt  HSM-99-72-30;  N°Vember'  1973‘  (also  supported  by  NIOSH 

a??SprhE‘  NaHI/  Hab,bi/  "Nonlinear  Adaptive  Recursive  Im- 
hanCemeno#  Proceedings  of  First  International  Joint 
PP?  485-491?"  PattGrn  Rec°nnition,  Washington,  D.C.,  Oct. 1973, 

p!ylllamiK*  ?raFt/  "D,Klta1  Transforms,"  1973  IEEE  NEREM  Sicnal 
Processing  Seminar,  Boston,  Mass.,  November,  1973.  g 

tiii'^Seve^h^A  ,7ra?Sf0rm  'TaRe  Cod,nG  Spectrum  Extrapola- 
ences,  January, m"  1  n ter na C 1  onal  Conference  on  System  Sci- 

t  ion1  EEEPTran^arM  rrel  3t  1  a*1  TecbnIdues  of  Image  Registra- 
May,  1974?E  Transactlons  on  Aerospace  and  Electronic  Systems, 


-152- 


1 


William  K.  Pratt,  Wen  H.  Chen,  Lloyd  R.  Welch, "Slant  Transform 
lmap,p  Codng,  IEEE  Transactions  on  Communi  cat  t  ons,  ( to  appear 


Alexander  A.  Sawchuk,  "Space- Va r i ant  Image  Restoration  by  Co¬ 
ordinate  Transformations",  Journal  of  the  Optical  Society  of 
America,  Vol.  64,  February,  1974,  pp.  124-130. 


James  Winter,  Mark  A.  Stein,  and  Alexander  A. 
verse  Axial  Tomography  Using  Coherent  Optical 

_  r  r  _ 


Sawchuk,  "Trans" 
verse  Axial  Tomography  Us  ins  Coherent  Optical  Image  Recon¬ 
struction,  Proceedings  of  Conference  on  Computerized  Trans¬ 
verse  Axial  Tomography,  UCLA,  Los  Angeles,  November,  1973. 


